基於有限體積法和交錯網格的SIMPLE演算法推導及實現

Mestro發表於2024-05-02

基於有限體積法和交錯網格的SIMPLE演算法推導及實現

SIMPLE演算法，半隱式速度壓力耦合演算法，是專門求解不可壓流體流動的演算法。由於不可壓流體控制方程中，密度常常被視為常數，沒有表徵流體密度、壓力、溫度聯絡的狀態方程，壓力以梯度項的形式存在於動量方程中，無法顯性表達或者直接求解，造成了求解上的困難。因此需要對速度和壓力進行解耦，其基本思想是，透過預設的壓力場，代入到動量方程中，求解各方向上的速度場；然後根據質量守恆方程構建壓力泊松方程，進行壓力修正，再進行速度修正，反覆迭代，直到達到預設的收斂條件。

交錯網格

與交錯網格相對的是同位網格，再RC插值方法提出之前，壓力、速度都存於同一套網格中，共用網格節點。例如x方向的壓力梯度項採用中心差分格式:

\[\frac{\partial P}{\partial x} = \frac{P_E - P_W}{2 \delta x} \]

並沒有使用到中心節點\(P_C\)的值，而在棋盤振盪的壓力場時，很容易被認為是無壓力梯度，即均勻壓力場。造成脫離真實物理的情況。因此交錯網格被提出，交錯網格中，壓力、溫度、密度這一類標量存在一套網格中，而各方向上的速度則存在另外一套網格中，速度網格的中心與壓力網格中心並不重合，存在一定偏差，速度網格中心一般定義在壓力網格面的中心上。如果是均勻網格，則速度網格與壓力網格偏差了“半個網格”。交錯網格的優勢很多，例如速度網格中心在壓力網格面上，在求解壓力修正方程時，不需要對速度再進行插值。避免了離散壓力梯度存在的“棋盤格”問題。但是交錯網格的致命缺陷就是使用非常複雜和繁瑣，儲存成本高，程式編制難度高，邊界條件比較難處理，都是致命缺陷。

穩態二維各向同性流體動力學方程

二維情況下穩態動力學方程組如下：

\[\frac\partial{\partial x}(\rho uu)+\frac\partial{\partial y}(\rho vu)=\frac\partial{\partial x}\Bigg(\mu\frac{\partial u}{\partial x}\Bigg)+\frac\partial{\partial y}\Bigg(\mu\frac{\partial u}{\partial y}\Bigg)-\frac{\partial p}{\partial x}+S_u \]

\[\frac\partial{\partial x}(\rho uv)+\frac\partial{\partial y}(\rho vv)=\frac\partial{\partial x}\Bigg(\mu\frac{\partial v}{\partial x}\Bigg)+\frac\partial{\partial y}\Bigg(\mu\frac{\partial v}{\partial y}\Bigg)-\frac{\partial p}{\partial y}+S_v \]

連續性方程:

\[\frac\partial{\partial x}(\rho u)+\frac\partial{\partial y}(\rho v)=0 \]

其中，對於各向同性的不可壓流體而言，密度\(\rho\)和粘性係數\(\mu\)為常數，可以提出微分符號外，並進一步化簡。使用交錯網格，即標量、向量（分量）分別存在一套網格中，彼此錯位。但是，在列離散方程時，以變數為核心，以變數所處的網格為依據進行離散。例如，對x方向的動量方程進行離散，所求解的是速度的x分量\(u\)，就以\(u\)網格為核心進行離散，所使用的網格編號就是\(u\)網格的網格編號。

上圖是交錯網格的示意圖，對應的是流場內部節點。速度u分量的網格中心(u-cell)，剛好位於標量p(p-cell)網格的左右邊界上，速度v分量的網格中心位於標量p網格的上下邊界上。交錯網格因為網格編號比較複雜，容易搞混，實際上只要記住一點，就可以理清楚三套網格的位置差異，即以行號i，列號j的p-cell網格為基準，行號i，列號j的u-cell網格在p-cell網格的左邊界，行號i，列號j的v-cell網格在p-cell網格的右邊界。記住這個轉換關係就可以理清在程式編寫時，同一i,j對於不同網格的相對位置關係。

動量方程離散

對對流項，採用一階迎風格式，擴散項採用前向差分格式，注意，因為速度的方向是未知的，並不一定沿著座標正方向，所以需要寫成通用格式離散。首先忽略體積力源項，採用有限體積法對方程左右進行二維情況下的體積分，將方程進行半離散，可得：

\[\int {{{\left. {\rho uu} \right|}_e}dy} - \int {{{\left. {\rho uu} \right|}_w}dy} + \int {{{\left. {\rho vu} \right|}_n}dx} - \int {{{\left. {\rho vu} \right|}_s}dx} = \left( {\int {{{\left. {\mu \frac{{\partial u}}{{\partial x}}} \right|}_e}dy} - \int {{{\left. {\mu \frac{{\partial u}}{{\partial x}}} \right|}_w}dy} } \right) + \left( {\int {\mu {{\left. {\frac{{\partial u}}{{\partial y}}} \right|}_n}dx} - \int {\mu {{\left. {\frac{{\partial u}}{{\partial y}}} \right|}_s}dx} } \right) + \int { - \frac{{\partial p}}{{\partial x}}dxdy} \]

將對流項和擴散項進行合併，再進行離散，以從u-cell右側介面流入網格的通量為例：

\[\int {{{\left. {\rho uu} \right|}_e}dy} - \int {{{\left. {\mu \frac{{\partial u}}{{\partial x}}} \right|}_e}dy} = \left( {\left( {\frac{{{{\left| {\rho {u_{n - 1}}} \right|}_e} + {{\left. {\rho {u_{n - 1}}} \right|}_e}}}{2}} \right){u_P} - \left( {\frac{{{{\left| {\rho {u_{n - 1}}} \right|}_e} + {{\left. {\rho {u_{n - 1}}} \right|}_e}}}{2}} \right){u_E} - \frac{{{u_{E - }}{u_P}}}{{\Delta x}}} \right)\Delta y \]

上式中，\(\rho u_{n-1}\)為質量通量，由於原本的\(\rho u u\)為非線性項，直接迭代較為困難。所以這裡採用延遲修正，即使用上一步迭代計算得到的速度，來計算質量通量。當迭代殘差越小時，上一步迭代的結果越接近當前迭代步計算得到的結果，實現一種類似“追趕”效果，直到兩者誤差在可接受範圍內。\(u_P\)為當前u-cell網格中的u速度，\(u_E\)為右側u-cell中心的u速度。\(\int {{{\left. {\rho uu} \right|}_e}dy}\)採用一階迎風格式進行離散。而\(\int {{{\left. {\mu \frac{{\partial u}}{{\partial x}}} \right|}_e}dy}\)則使用了前向差分，一階格式進行離散得到。這裡\(\mu\)被認為是一個常數，簡化了運算，否則，\(\mu\)實際上需要被當做一個場量來處理，需要透過緊鄰點的\(\mu\)值得到。

注意\(\rho u_{n-1}\)是u-cell右側邊界上的通量值，它仍然需要透過插值來獲得，如下圖所示：

透過對介面鄰近網格的體心值進行插值可以得到介面上的值：

\[\rho u_{n-1} = \rho\frac{u_P^{n-1} + u_E^{n-1}}{2} \]

因為流體是均質不可壓流體，密度被設為常數。

同理可得其餘項的離散方程：

\[\int {{{\left. {\rho uu} \right|}_w}dy} - \int {{{\left. {\mu \frac{{\partial u}}{{\partial x}}} \right|}_w}dy} = \left( {\left( {\frac{{{{\left| {\rho {u_{n - 1}}} \right|}_w} + {{\left. {\rho {u_{n - 1}}} \right|}_w}}}{2}} \right){u_W} - \left( {\frac{{{{\left| {\rho {u_{n - 1}}} \right|}_e} - {{\left. {\rho {u_{n - 1}}} \right|}_e}}}{2}} \right){u_P} - \frac{{{u_{P - }}{u_W}}}{{\Delta x}}} \right)\Delta y \]

\[\int {{{\left. {\rho vu} \right|}_n}dy} - \int {{{\left. {\mu \frac{{\partial u}}{{\partial y}}} \right|}_n}dx} = \left( {\left( {\frac{{{{\left| {\rho {v_{n - 1}}} \right|}_n} + {{\left. {\rho {v_{n - 1}}} \right|}_n}}}{2}} \right){u_P} - \left( {\frac{{{{\left| {\rho {v_{n - 1}}} \right|}_n} - {{\left. {\rho {v_{n - 1}}} \right|}_n}}}{2}} \right){u_P} - \frac{{{u_{N - }}{u_P}}}{{\Delta y}}} \right)\Delta x \]

\[\int {{{\left. {\rho vu} \right|}_s}dy} - \int {{{\left. {\mu \frac{{\partial u}}{{\partial y}}} \right|}_s}dx} = \left( {\left( {\frac{{{{\left| {\rho {v_{n - 1}}} \right|}_s} + {{\left. {\rho {v_{n - 1}}} \right|}_s}}}{2}} \right){u_S} - \left( {\frac{{{{\left| {\rho {v_{n - 1}}} \right|}_s} - {{\left. {\rho {v_{n - 1}}} \right|}_s}}}{2}} \right){u_P} - \frac{{{u_{P - }}{u_S}}}{{\Delta y}}} \right)\Delta x \]

這裡，涉及到\(\rho v_{n-1}\)，需要找到對應的v-cell網格進行插值得到，如下圖所示：

藍色為v-cell，紅色為u-cell，u-cell的上下介面的質量通量，，假定\(u_P\)的網格索引號為i,j，則對應到v-cell網格中的索引如下，\(\rho v_{n-1}\)透過相鄰v-cell的體心值插值獲得:

\[\rho v_{n-1}|_{s}=\rho\frac{v_P+v_W}{2}=\rho\frac{v_{i,j}+v_{i,j-1}}{2} \]

\[\rho v_{n-1}|_{n}=\rho\frac{v_N+v}{2}=\rho\frac{v_{i+1,j}+v_{i+1,j-1}}{2} \]

最後是壓力梯度項的離散，這裡需要說明，雖然此處是對壓力梯度進行離散，但是求解的變數為u，以u-cell的位置為基準進行離散，這裡體現出一個交錯網格的優勢，它十分巧妙的讓幾套網格的中心落在邊界上。u-cell的中心，其實是壓力網格的左邊界。因此以u-cell為中心的壓力梯度離散，可以直接採用p-cell的體心值前向差分得到，不需要像上述過程在邊界處再插值：

\[\int \frac{\partial p}{\partial x}dxdy= (P_P-P_W)\Delta y \]

同理對y方向上的動量方程也進行上述離散過程，不再贅述。對上述格式進行整理可以獲得:

\[a_{{i,\tilde{f}}}u_{{i,\tilde{f}}}=\sum a_{nb}u_{nb}-({p_{I,\tilde{f}}-p_{I-1,\tilde{f}}})\Delta y \]

\[a_{I,j}v_{I,j}=\sum a_{nb}v_{nb}+(p_{I,j-1}-p_{I,j})\Delta x \]

SIMPLE演算法速度和壓力修正

接下來是SIMPLE演算法的核心步驟：

假定\(u^*\)和\(v^*\)為某一迭代步計算後得到的解，則有(\(A_{ij}\)即\(\Delta x\)、\(\Delta y\))：

\[a_{i,f}u_{i,f}^*=\sum a_{nb}u_{nb}^*+(p_{I-1,f}^*-p_{I,f}^*)A_{i,f}+b_{i,f} \]

\[a_{I,j}v_{I,j}^*=\sum a_{nb}v_{nb}^*+(p_{I,j-1}^*-p_{I,j}^*)A_{I,j}+b_{I,j} \]

將上式與假定的收斂解做差，可得:

\[a_{{i,j}}(u_{{i,j}}-u_{{i,j}}^{*})=\sum a_{nb}(u_{nb}-u_{nb}^{*})+[(p_{{I-1,j}}-p_{{I-1,j}}^{*})-(p_{{I,j}}-p_{{I,j}}^{*})]A_{{i,j}} \]

\[a_{I,j}(v_{I,j}-v_{I,j}^{*})=\sum a_{nb}(v_{nb}-v_{nb}^{*})+[(p_{I,j-1}-p_{I,j-1}^{*})-(p_{I,j}-p_{I,j}^{*})]A_{I,j} \]

進一步簡寫可以得到:

\[a_{i,j}u_{i,j}^{\prime}=\sum a_{nb}u_{nb}^{\prime}+(p_{I-1,j}^{\prime}-p_{I,j}^{\prime})A_{i,j}\\a_{I,j}v_{I,j}^{\prime}=\sum a_{nb}v_{nb}^{\prime}+(p_{I,j-1}^{\prime}-p_{I,j}^{\prime})A_{I,j} \]

上式中，\(u_{i,j}^{\prime}\)和\(v_{I,j}^{\prime}\)則分別是u-cell和v-cell速度修正值，因此速度的更新，則由迭代值\(u^*\)和\(u^{\prime}\)決定：

\[u=u^*+u^{\prime} \]

\[v=v^{*}+v^{\prime} \]

上式中，\(\sum a_{nb}(u_{nb}-u_{nb}^{*})\)為相鄰網格節點對速度修正值的影響。SIMPLE演算法在迭代過程中，將這一值忽略，即速度修正僅與壓力修正值有關。這樣做，近鄰點無法顯性表達，直接迭代處理較為麻煩，二是，隨著速度逐漸收斂，近鄰點的影響逐漸減弱，直到“忽略不計”，並不會影響最終的結果：

\[u_{i,j}^{\prime}=d_{i,j}(p_{I-1,j}^{\prime}-p_{I,j}^{\prime})\\v_{I,j}^{\prime}=d_{I,j}(p_{I,j-1}^{\prime}-p_{I,j}^{\prime}) \]

\[u_{i,j}=u_{i,j}^*+d_{i,j}(p_{i-1,j}^\prime-p_{i,j}^\prime) \]

\[v_{I,j}=v_{I,j}^*+d_{I,j}(p_{I,j-1}^{\prime}-p_{I,j}^{\prime}) \]

到這一步，我們終於把壓力修正值與速度修正值聯絡起來了。可是僅靠動量方程，我們還是無法得到將兩者求解出來。這時，還剩下連續性方程，連續性方程沒有壓力錶達式的參與，僅由速度和密度決定。更像是一個準確性判據，即速度必須需要滿足連續性方程，才是符合物理解的方程。

\[[(\rho uA)_{i+1,j}-(\rho uA)_{i,j}]+[(\rho vA)_{I,j+1}-(\rho vA)_{I,j}]=0 \]

將u和v的迭代表示式代入上述方程，可以得到壓力泊松方程：

\[\begin{aligned}&[\rho_{i+1,j}A_{i+1,j}(u_{i+1,j}^{*}+d_{i+1,j}(p_{I,j}^{\prime}-p_{I+1,j}^{\prime}))-\rho_{i,j}A_{i,j}(u_{i,j}^{*}\\&+d_{i,j}(p_{I-1,j}^{\prime}-p_{I,j}^{\prime}))]+[\rho_{I,j+1}A_{I,j+1}(v_{I,j+1}^{*}+d_{I,j+1}(p_{I,j}^{\prime}-p_{I,j+1}^{\prime}))\\&-\rho_{I,j}A_{I,j}(v_{I,j}^{*}+d_{I,j}(p_{I,j-1}^{\prime}-p_{I,j}^{\prime}))]=0\end{aligned} \]

進一步化簡得到：

\[a_{I,j}p_{I,j}^{\prime}=a_{I+1,j}p_{I+1,j}^{\prime}+a_{I-1,j}p_{I-1,j}^{\prime}+a_{I,j+1}p_{I,j+1}^{\prime}+a_{I,j-1}p_{I,j-1}^{\prime}+b_{I,j}^{\prime} \]

壓力泊松方程，求解的是各個節點上的壓力修正值。仍然是需要迭代求解，並非顯性可以直接遞推獲得全域的壓力修正值。

邊界條件

邊界條件應該是交錯網格中，最容易讓人誤判的地方，參考書目中對邊界條件講的並不太多。首先，這是與交錯網格劃分方式相關。一般而言，對於真實流體域劃分網格，以標量網格為核心，速度網格再錯開。如下圖所示：

以壓力網格為核心，速度網格u和v相對於p網格中心進行偏移。對於無滑移邊界條件而言，邊界上的速度值已經給定，而u和v的邊界網格中心，剛好落在物理邊界上。所以直接對速度網格賦予邊界值，在求解壓力泊松方程。

而對於p網格，其邊界條件的確定，則對應到壓力泊松方程，具體參考汪洋博士的開源文件（B站名：大官人學CFD）。

本文所使用的網格劃分方式，與上述方式略有不用，更"類似"於有限差分法的網格劃分方式，即p-cell的中心在物理邊界上，導致邊界上的速度網格向外延伸，超出邊界成為虛擬網格。

對於u網格（上圖中紅色網格）：

u網格在左右邊界向外延伸出一個虛擬網格節點，而在上下邊界中，則u-cell中心落在邊界上，因此對於無滑移邊界或者速度給定的邊界條件，則，上下邊界的u-cell網格可以直接給定速度值。而對於左右邊界，超出實際物理邊界的虛擬邊界網格的賦值，則認為與內層節點互為相反數（對於無滑移邊界）。因為物理邊界實際處於虛擬節點與內層節點的中心位置，u=0，所以根據插值關係得到虛擬節點的賦值與內層節點的計算值互為相反數。

對於v網格（上圖中藍色網格）：

與u網格處理方式極為相似，不同的是，計算域左右物理邊界，剛好處於v網格的中心位置，所以可以直接賦值。而上下邊界處，v網格向外延伸出一層虛擬節點，和u網格處理方法一致，在無滑移邊界條件下，邊界虛擬網格節點的值與內層網格節點值互為相反數。

對於p網格（上圖中黑色網格）：

對於p網格，要分情況討論，還是要落到原本的壓力修正方程中。四個角點的p網格，首先考慮到，壓力是一個相對意義大於絕對意義的物理量，CFD計算中關注壓差，而非絕對壓力，壓差是驅動流體運動的關鍵。計算壓力需要設定一個壓力參考點，由於上邊界設定為速度邊界，所以習慣性設定左下角角點處的壓力為0，作為一個計算參考值。

緊接著對於其餘3個角點，將網格邊界上的通量代入到原本壓力修正方程中，可以得到，對於右下角角點，AE=0，AS=0。左上角角點，AW=0，AN=0。右上角角點，AN=0，AE=0。

然後是除角點外的邊界網格，對於左邊界，AW=0；右邊界，AE=0；上邊界，AN=0；下邊界, AS = 0。

最後是SIMPLE演算法的流程：

計算程式

以二維驅動方腔流動為例，編寫程式，實現使用有限體積法和交錯網格計算方腔內的速度和壓力分佈（個人覺得最難的部分，其實還是在於邊界條件的理解上）：

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Apr 28 09:46:16 2024

@author: mestro
"""
from math import *
import numpy as np
import matplotlib.pyplot as plt


L = 1.0
N = 51
dx = L/(N-1)
dy = L/(N-1)
velocity = 2.0
rho = 1
miu = 0.01

P = np.zeros([N,N])
U = np.zeros([N,N+1])
V = np.zeros([N+1,N])
U[N-1,:] = velocity
dU = np.zeros([N,N+1])
dV = np.zeros([N+1,N])


P_iter = P.copy()
U_iter = U.copy()
V_iter = V.copy()

NN = N*N
iter = 0
err = 1.0

while (iter < 1000) and (err > 1e-5):
    for i in range(1,N-1):
        for j in range(1,N):
            rho_u_e = 0.5*(U[i,j] + U[i,j+1])*rho;
            rho_u_w = 0.5*(U[i,j] + U[i,j-1])*rho;
            rho_v_s = 0.5*(V[i,j] + V[i,j-1])*rho;
            rho_v_n = 0.5*(V[i+1,j] + V[i+1,j-1])*rho;
            
            AE = 0.5*(abs(rho_u_e) + rho_u_e)*dy + miu*dy/dx;
            AW = 0.5*(abs(rho_u_w) - rho_u_w)*dy + miu*dy/dx;
            AN = 0.5*(abs(rho_v_n) + rho_v_n)*dx + miu*dx/dy;
            AS = 0.5*(abs(rho_v_s) - rho_v_s)*dx + miu*dx/dy;
            
            AEE = 0.5*(abs(rho_u_e) - rho_u_e)*dy + miu*dy/dx;
            AWW = 0.5*(abs(rho_u_w) + rho_u_w)*dy + miu*dy/dx;
            ANN = 0.5*(abs(rho_v_n) - rho_v_n)*dx + miu*dx/dy;
            ASS = 0.5*(abs(rho_v_s) + rho_v_s)*dx + miu*dx/dy;
            
            Ap = (AE+AW+AN+AS);
            U_iter[i,j] = 1/Ap*(AEE*U[i,j+1] + AWW*U[i,j-1] + ANN*U[i+1,j] + ASS*(U[i-1,j]) - (P[i,j] - P[i,j-1])*dy);
            dU[i,j] = dy/Ap;
    #bottom mesh
    i = 0;
    for j in range(1,N):
        rho_u_e = 0.5*(U[i,j] + U[i,j+1])*rho;
        rho_u_w = 0.5*(U[i,j] + U[i,j-1])*rho;
        #rho_v_s = 0.5*(V[i,j] + V[i,j-1])*rho;
        rho_v_n = 0.5*(V[i+1,j] + V[i+1,j-1])*rho;
        
        AE = 0.5*(abs(rho_u_e) + rho_u_e)*dy + miu*dy/dx;
        AW = 0.5*(abs(rho_u_w) - rho_u_w)*dy + miu*dy/dx;
        AN = 0.5*(abs(rho_v_n) + rho_v_n)*dx + miu*dx/dy;
        #AS = 0.5*(abs(rho_v_s) - rho_v_s)*dx + miu*dx/dy;
        As = 0;
        Ap = (AE+AW+AN+AS);
        dU[i,j] = dy/Ap;
    #top mesh
    i = N-1
    for j in range(1,N):
        rho_u_e = 0.5*(U[i,j] + U[i,j+1])*rho;
        rho_u_w = 0.5*(U[i,j] + U[i,j-1])*rho;
        rho_v_s = 0.5*(V[i,j] + V[i,j-1])*rho;
        #rho_v_n = 0.5*(V[i+1,j] + V[i+1,j-1])*rho;
        
        AE = 0.5*(abs(rho_u_e) + rho_u_e)*dy + miu*dy/dx;
        AW = 0.5*(abs(rho_u_w) - rho_u_w)*dy + miu*dy/dx;
        #AN = 0.5*(abs(rho_v_n) + rho_v_n)*dx + miu*dx/dy;
        AN = 0.0;
        AS = 0.5*(abs(rho_v_s) - rho_v_s)*dx + miu*dx/dy;
        Ap = (AE+AW+AN+AS);
        dU[i,j] = dy/Ap;
    #Apple BCs
    U_iter[:,0] = -U_iter[:,1]; #left
    U_iter[1:N-1,N] = -U_iter[1:N-1,N-1]; #right
    U_iter[0,:] = 0.0 #bottom
    U_iter[N-1,:] = velocity #top
    
    
    # V equation
    for i in range(1,N):
        for j in range(1,N-1):
            rho_u_e = 0.5*(U[i,j] + U[i-1,j])*rho;
            rho_u_w = 0.5*(U[i,j+1] + U[i-1,j+1])*rho;
            rho_v_n = 0.5*(V[i,j] + V[i+1,j])*rho;
            rho_v_s = 0.5*(V[i,j] + V[i-1,j])*rho;
            
            AE = 0.5*(abs(rho_u_e) + rho_u_e)*dy + miu*dy/dx;
            AW = 0.5*(abs(rho_u_w) - rho_u_w)*dy + miu*dy/dx;
            AN = 0.5*(abs(rho_v_n) + rho_v_n)*dx + miu*dx/dy;
            AS = 0.5*(abs(rho_v_s) - rho_v_s)*dx + miu*dx/dy;
            
            AEE = 0.5*(abs(rho_u_e) - rho_u_e)*dy + miu*dy/dx;
            AWW = 0.5*(abs(rho_u_w) + rho_u_w)*dy + miu*dy/dx;
            ANN = 0.5*(abs(rho_v_n) - rho_v_n)*dx + miu*dx/dy;
            ASS = 0.5*(abs(rho_v_s) + rho_v_s)*dx + miu*dx/dy;
            
            Ap = (AE+AW+AN+AS);
            V_iter[i,j] = 1.0/Ap*(AEE*V[i,j+1] + AWW*V[i,j-1] + ANN*V[i+1,j] + ASS*V[i-1,j] - (P[i,j] - P[i-1,j])*dx);
            dV[i,j] = dx/Ap;
    #left
    j = 0
    for i in range(1,N):
        rho_u_e = 0.5*(U[i,j] + U[i-1,j])*rho;
        #rho_u_w = 0.5*(U[i,j+1] + U[i-1,j+1])*rho;
        rho_v_n = 0.5*(V[i,j] + V[i+1,j])*rho;
        rho_v_s = 0.5*(V[i,j] + V[i-1,j])*rho;
        
        AE = 0.5*(abs(rho_u_e) + rho_u_e)*dy + miu*dy/dx;
        #AW = 0.5*(abs(rho_u_w) - rho_u_w)*dy + miu*dy/dx;
        AW = 0.0
        AN = 0.5*(abs(rho_v_n) + rho_v_n)*dx + miu*dx/dy;
        AS = 0.5*(abs(rho_v_s) - rho_v_s)*dx + miu*dx/dy;
        Ap = (AE+AW+AN+AS);
        dV[i,j] = dx/Ap;
    
    #right0, L, N
    j = N-1
    for i in range(1,N):
        #rho_u_e = 0.5*(U[i,j] + U[i-1,j])*rho;
        rho_u_w = 0.5*(U[i,j+1] + U[i-1,j+1])*rho;
        rho_v_n = 0.5*(V[i,j] + V[i+1,j])*rho;
        rho_v_s = 0.5*(V[i,j] + V[i-1,j])*rho;
        
        #AE = 0.5*(abs(rho_u_e) + rho_u_e)*dy + miu*dy/dx;
        AE = 0.0;
        AW = 0.5*(abs(rho_u_w) - rho_u_w)*dy + miu*dy/dx;
        AN = 0.5*(abs(rho_v_n) + rho_v_n)*dx + miu*dx/dy;
        AS = 0.5*(abs(rho_v_s) - rho_v_s)*dx + miu*dx/dy;
        Ap = (AE+AW+AN+AS);
        dV[i,j] = dx/Ap;
    #Apply BCs

    V_iter[:,0] = 0.0;
    V_iter[:,N-1] = 0.0;
    V_iter[0,1:N-1] = -V_iter[1,1:N-1];
    V_iter[N,1:N-1] = -V_iter[N-1,1:N-1];
    
    U_old = U.copy();
    V_old = V.copy();
    bp = np.zeros([NN,1]);
    #pressure fix
    for i in range(N):
        for j in range(N):
            index = i*N+j;
            bp[index] = (rho*U_iter[i,j]*dy - rho*U_iter[i,j+1]*dy + rho*V_iter[i,j]*dx - rho*V_iter[i+1,j]*dx);
    bp[0] = 0.0;
    
    APP = np.zeros([NN,NN]);
    #left bottom
    i = 0
    j = 0
    index = i*N + j
    # Ae = -rho*dU[i,j+1]*dy;
    # An = -rho*dV[i+1,j]*dx;
    # Ap = -(Ae + An);
    # APP[index,index+1] = Ae;
    # APP[index,index+N] = An;
    # APP[index,index] = Ap;
    APP[index,index] = 1;
    
    #right bottom
    i = 0;
    j = N-1;
    index = i*N + j
    Aw = -rho*dU[i,j]*dy;
    An = -rho*dV[i+1,j]*dx;
    Ap = -(Aw + An);
    APP[index,index - 1] = Aw;
    APP[index,index + N] = An;
    APP[index,index] = Ap;
    
    #left top
    i = N-1
    j = 0;
    index = i*N + j
    As = -rho*dV[i,j]*dx;
    Ae = -rho*dU[i,j+1]*dy;
    Ap = -(As + Ae);
    APP[index,index+1] = Ae;
    APP[index,index-N] = As;
    APP[index,index] = Ap;
    
    #right top
    i = N-1
    j = N-1
    index = i*N+j
    Aw = -rho*dU[i,j]*dy;
    As = -rho*dV[i,j]*dx;
    Ap = -(Aw + As);
    APP[index,index] = Ap
    APP[index,index-1] = Aw
    APP[index,index-N] = As 

    i = 0;
    for j in range(1,N-1):
        index = i*N+j;
        Aw = -rho*dU[i,j]*dy;
        An = -rho*dV[i+1,j]*dx;
        Ae = -rho*dU[i,j+1]*dy;
        Ap = -(Aw + An + Ae);
        APP[index,index] = Ap;
        APP[index,index-1] = Aw;
        APP[index,index+N] = An;
        APP[index,index+1] = Ae;
        
        
        
    i = N-1;
    for j in range(1,N-1):
        index = i*N+j
        Aw = -rho*dU[i,j]*dy;
        As = -rho*dV[i,j]*dx;
        Ae = -rho*dU[i,j+1]*dy;
        Ap = -(Aw + As + Ae);
        APP[index,index] = Ap;
        APP[index,index - 1] = Aw;
        APP[index,index - N] = As;
        APP[index,index + 1] = Ae;
    
    j = 0;
    for i in range(1,N-1):
        index = i*N+j;
        Ae = -rho*dU[i,j+1]*dy;
        An = -rho*dV[i+1,j]*dx;
        As = -rho*dV[i,j]*dx;
        Ap = -(Ae + An + As);
        APP[index,index] = Ap;
        APP[index,index + 1] = Ae;
        APP[index,index + N] = An;
        APP[index,index - N] = As;
    
    j = N-1;
    for i in range(1,N-1):
        index = i*N + j;
        Aw = -rho*dU[i,j]*dy;
        An = -rho*dV[i+1,j]*dx;
        As = -rho*dV[i,j]*dx;
        Ap = -(Aw + An + As);
        APP[index,index] = Ap;
        APP[index,index - 1] = Aw;
        APP[index,index + N] = An;
        APP[index,index - N] = As;
    
    for i in range(1,N-1):
        for j in range(1,N-1):
            index = i*N + j
            Aw = -rho*dU[i,j]*dy;
            An = -rho*dV[i+1,j]*dx;
            As = -rho*dV[i,j]*dx;
            Ae = -rho*dU[i,j+1]*dy;
            Ap = -(Aw + An + As + Ae);
            APP[index,index] = Ap;
            APP[index,index - 1] = Aw;
            APP[index,index + N] = An;
            APP[index,index - N] = As;
            APP[index,index + 1] = Ae;
            
            
    # pressure correction
    p_fix = np.linalg.solve(APP, bp)
    P_fix_matrix = np.zeros([N,N]);
    for i in range(N):
        for j in range(N):
            index = i*N+j
            P[i,j] = 0.3*p_fix[index] + P[i,j];
            P_fix_matrix[i,j] = p_fix[index];
    P[0,0] = 0.0;
    #velocity update
    for i in range(1,N-1):
        for j in range(1,N):
            U[i,j] = U_iter[i,j] + dU[i,j]*(P_fix_matrix[i,j-1] - P_fix_matrix[i,j]);
    for i in range(1,N):
        for j in range(1,N-1):
            V[i,j] = V_iter[i,j] + dV[i,j]*(P_fix_matrix[i-1,j] - P_fix_matrix[i,j]);
    
    
    #Apple BCs
    U[1:N-1,0] = -U[1:N-1,1]; #left
    U[1:N-1,N] = -U[1:N-1,N-1]; #right
    U[0,:] = 0.0 #bottom
    U[N-1,:] = velocity #top
    
    V[0,1:N-1] = -V[1,1:N-1];
    V[N,1:N-1] = -V[N-1,1:N-1];
    V[:,0] = 0.0;
    V[:,N-1] = 0.0;
    
    err1 = np.max(np.abs(U-U_old))
    err2 = np.max(np.abs(V-V_old))
    err = max(err1,err2)
    
    iter = iter + 1
    print("the iter num is " + str(iter) + " and max err is " + str(err) + "\n")
    


#plot
x = np.linspace(dx/2, 1 - dx/2,N-1);
y = np.linspace(0, 1, N)
X, Y = np.meshgrid(x, y)  # Generate 2D grid coordinates
plt.figure()
plt.contourf(X,Y, U[:,1:N], levels=20, cmap='jet')  # Adjust levels and cmap as needed
plt.colorbar(label='Velocity UX')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Velocity Contour Plot')
plt.show()


x = np.linspace(0,1,N);
y = np.linspace(0, 1, N)
X, Y = np.meshgrid(x, y)  # Generate 2D grid coordinates
plt.figure()
plt.contourf(X,Y,P, levels=20, cmap='jet')  # Adjust levels and cmap as needed
plt.colorbar(label='pressure P')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Pressure')
plt.show()

小結

本文主要是使用Python實現了基於有限體積法和交錯網格的SIMPLE演算法，對流項採用一階迎風格式，擴散項採用前向差分，其實離散過程並不複雜，最複雜和最難分清楚的還是在於邊界條件，個人覺得這也是書中講述比較少的地方。文中有錯誤非常歡迎指出，尤其是理解上和原則性錯誤。文章會同步到個人部落格（https://bihengx.github.io）

參考資料：

B站曾導SJTU的CFD從0到1的課程
B站大官人學CFD，汪洋博士的開源OpenFOAM基礎課程講義
參考書：An Introduction to Computational Fluid Dynamics The Finite Volume Method 2nd Edition.

CRC演算法原理、推導及實現
2024-08-16
演算法
模擬積體電路設計系列部落格——7.5.6 時間交錯型ADC
2024-06-27
Etsy使用交錯新演算法實現更快的ML實驗
2022-08-15
演算法
simple-mybatis的原理及實現
2018-06-19
MyBatis
數值積分公式及龍貝格（Romberg）演算法實現matlab
2020-10-13
公式演算法Matlab
有限元模擬有限體積模擬
2024-04-05
基於卷積神經網路和tensorflow實現的人臉識別
2018-04-01
卷積神經網路
深度學習基礎-基於Numpy的卷積神經網路（CNN）實現
2022-06-09
深度學習卷積神經網路CNN
基於MySql和Sails.js的RESTful風格的api實現
2018-07-03
MySqlAIJSRESTAPI
【Pytorch】基於卷積神經網路實現的面部表情識別
2021-12-17
PyTorch卷積神經網路
Git基礎及Simple Gitflow
2021-08-04
Git
我是如何基於二階段遞交及悲觀鎖實現分散式事務的
2019-07-14
分散式
推導有限元格式的幾種方法6
2024-10-07
關於個人部落格的選擇及實現
2022-01-17
基於libcurl實現REST風格http/https的get和post
2018-10-16
RESTHTTP
卷積神經網路的原理及Python實現
2024-05-25
卷積神經網路Python
卷積神經網路(CNN)反向傳播演算法公式詳細推導
2020-04-04
卷積神經網路CNN反向傳播演算法公式
【基於物品的推薦演算法】
2018-01-28
演算法
Survey | 基於圖卷積網路的藥物發現方法
2019-09-04
卷積
貝塞爾曲線原理、推導及Matlab實現
2024-07-04
Matlab
Go實戰 | 基於本地記憶體的快取的應用及實現
2021-11-20
Go記憶體快取
基於Redis實現線上遊戲積分排行榜
2020-09-23
Redis遊戲
基於BKDRhash實現Hash演算法
2014-09-12
演算法
TF2.keras 實現基於卷積神經網路的影象分類模型
2020-03-11
TF2Keras卷積神經網路模型
TF2.keras 實現基於卷積神經網路的影像分類模型
2020-03-11
TF2Keras卷積神經網路模型
靜態網格體和骨架網格體的區別
2024-04-06
經典論文復現 | 基於深度卷積網路的影象超解析度演算法
2018-12-14
卷積演算法
經典論文復現 | 基於深度卷積網路的影像超解析度演算法
2018-12-14
卷積演算法
[詳細推導]基於EKF的小車運動模型的python程式設計實現
2020-11-04
模型Python程式設計
Python 列表推導及優先順序佇列的實現
2019-01-10
Python佇列
基於Python和Tensorflow的電影推薦演算法
2018-05-16
Python演算法
TF2.keras 實現基於深度可分離卷積網路的影象分類模型
2020-03-11
TF2Keras卷積模型
TF2.keras 實現基於深度可分離卷積網路的影像分類模型
2020-03-11
TF2Keras卷積模型
基於軟體實現網站負載均衡(1) (轉)
2007-12-26
網站負載
演算法導論C語言實現: 演算法基礎
2013-06-30
演算法C語言
實現基於內容的電影推薦系統—程式碼實現
2024-04-07
最全面的圖卷積網路GCN的理解和詳細推導，都在這裡了!
2021-02-02
卷積GC
基於android的智慧導診的設計與實現
2019-02-12
Android

基於有限體積法和交錯網格的SIMPLE演算法推導及實現

基於有限體積法和交錯網格的SIMPLE演算法推導及實現

交錯網格

穩態二維各向同性流體動力學方程

動量方程離散

SIMPLE演算法速度和壓力修正

邊界條件

計算程式

小結

相關文章