導語

 

Paritosh 是 Wolfram 的核心開發(fā)人員，利用業(yè)余時間使用 Mathematica 來研究并模擬流體動力學問題，開發(fā)了WindTunnel2DLBM 程序包(https://blog.wolfram.com/data/uploads/2019/10/WindTunnel2DLBM.zip) 。LBM 與 IBM 的結(jié)合使用，對研究和分析流體流動是一個很好的工具。借助 Mathematica 的內(nèi)置函數(shù)，實現(xiàn)數(shù)字風洞的組裝變得非常簡單。

 

我在學生時代學習流體動力學時，學的都是復雜的方程式以及各種簡化和處理這些方程式的方法，以獲得某個結(jié)果。遺憾的是，當要獲得流體在不同情況下如何流動的直觀感受時，這讓想象力幾乎沒有什么空間發(fā)揮。當我上完第一節(jié)實驗流體動力學課后，我了解了如何使用不同的可視化技術(shù)來定性地了解流體行為。這些可視化方式提供了一種創(chuàng)造性地查看流體的方法，并且，效果令人驚艷。所有這些實驗和可視化都在風洞 (wind tunnel) 內(nèi)進行的。

創(chuàng)建我們的計算型風洞

 

風洞是實驗流體力學研究人員用來研究物體在空氣（或任何其他流體）中移動時所產(chǎn)生的影響的設(shè)備。由于物體本身無法在管道內(nèi)移動，因此通常采用高速風扇提供可控氣流，而將物體放置在氣流路徑中，從而產(chǎn)生與物體在靜止空氣中移動相同的效果。這種實驗對于理解物體的空氣動力學非常有用。

 

風洞有不同的類型。最簡單的風洞是空心管或方盒，該風洞的一端裝有風扇，另一端是敞開的，這種風洞稱為開放式回風風洞。從實驗上講，這種風洞不是最有效或最可靠的，卻很適合構(gòu)建一個計算型風洞（這正是我們這篇博文要實現(xiàn)的目標）。這是我們要開發(fā)的風洞的基本示意圖：

我們的風洞是一個二維風洞，流體從左側(cè)進入，然后從右側(cè)流出，頂部和底部是堅固的墻壁。應(yīng)該指出，由于這是一個計算風洞，我們在選擇邊界類型上有很大的靈活性。例如，我們可讓左壁、右壁和底壁固定，而頂壁可以移動，這樣我們就擁有了一個很常見的風洞類型。

 

當人們思考計算流體力學時，總會立刻想到著名的納維-斯托克斯方程。這組方程是決定流體流動行為的控制方程。在這一點上，我們似乎要使用納維-斯托克斯方程來幫助我們建造一個風洞。但事實證明，還有其他方法可以研究流體流動的行為，而無需求解納維-斯托克斯方程。格子玻爾茲曼方法 （LBM）就是這些方法中的一個。

 

如果使用納維-斯托克斯方程，我們面對的將是一個復雜的偏微分方程組（PDE）。為了在數(shù)值上求解，我們必須應(yīng)用各種技巧來離散化導數(shù)。一旦離散化完成，得到的將是一個龐大的非線性代數(shù)方程組有待求解。好累！而使用 LBM 方法，將完全繞過傳統(tǒng)方法。在 LBM 中不需要求解方程組。此外，許多運算（我稍后將介紹）完全是本地運算。這使得 LBM 成為高度并行的方法。

 

我們可以在一個非常簡單的框架中，將LBM視為一種"自下而上"的方法。在這種方法中，模擬在微觀或"格子"域中進行。假設(shè)您有一個物理域或宏觀域。我們要放大這個宏觀域中的單個點，將有一些粒子根據(jù)粒子之間相互作用的"規(guī)則"而相互作用：

例如，如果兩個粒子相互撞擊，它們?nèi)绾畏磸棧窟@些粒子遵循某種離散規(guī)則?，F(xiàn)在，如果我們根據(jù)這些規(guī)則讓粒子隨時間演變（可以看到這與元胞自動機密切相關(guān)），并獲取平均值，那么這些平均值可以用來描述某些宏觀量。例如，我們在上圖中看到的 HPP 模型（以哈代、波莫和德帕齊斯命名）可用于模擬氣體擴散。

 

盡管這種離散方法聽起來不錯（20世紀70年代中期的研究人員確實嘗試了其可行性），但它有許多缺點。其中一個主要問題是最終結(jié)果中有統(tǒng)計噪聲。正是由于這些原則問題和試圖解決這些問題的努力，LBM方法出現(xiàn)了。關(guān)于該方法的理論方面，網(wǎng)上有很多推導至最終方程的鏈接。在本文中，我想關(guān)注的是玻爾茲曼模擬的最終基礎(chǔ)機制，而不是理論方面。因此，我僅談一下開發(fā)風洞所需的最后方程。首先，假設(shè)密度和速度由以下方程描述：

……其中 fi 為分布函數(shù)，ex, i  和 ey, i 為離散速度，且具有下列值：

            
密度和速度的計算可以如下所示：

這九個離散速度與其各自的分布函數(shù) fi 相關(guān)聯(lián)，可以形象地用下圖表示：

在格子域中的每個（離散）點上，將存在九個分布函數(shù)。使用這九個離散速度和函數(shù) fi 的模型稱為 D2Q9 模型。如果分布函數(shù) fi 已知，則速度場已知。由于速度場隨空間和時間演化，我們認為這些分布函數(shù)也會隨空間和時間演化。控制 fi 的時空演化方程由下式給出：

這里的 Ωi 是一個復雜的"碰撞"項，它基本上決定了各種 fi 如何相互交互。經(jīng)過各種簡化和近似，此方程可以化簡為下式：

……其中 fieq 稱為平衡分布函數(shù)，τ 稱為松弛參數(shù)，δtLBM = 1 是格子玻爾茲曼域中的時間步長。該近似稱為 BGK 近似，生成的模型稱為格子 BGK 模型。關(guān)于這兩個項的詳細定義將在稍后介紹。這些函數(shù)的空間和時間演化分兩個步驟完成：（a） 遷移步驟；和 （b） 碰撞步驟。

 

由于δtLBM = 1 且

均為1或者0，遷移步驟由下式給出：

要將遷移步驟可視化，想象一下離散化的域。此域的每個網(wǎng)格點上有 9 個分布函數(shù)（如下圖所示）。請注意每個網(wǎng)格點的各種顏色。每個箭頭的長度指示相應(yīng)的分布函數(shù)的大小：

根據(jù)數(shù)學公式，每個分布將按以下方向進行遷移：

注意顏色和它們結(jié)束的位置。將注意力集中在中心點會有幫助。在遷移步驟之前，箭頭（表示不同的fi）都是綠色的。在遷移步驟之后，請注意綠色箭頭落在周圍網(wǎng)格點的位置。簡言之，這就是遷移步驟。在 Mathematica 中，此步驟使用內(nèi)置函數(shù)RotateLeft 或 RotateRight 可以非常容易地完成。這里我們使用 RotateRight：

當進行遷移時，必須特別注意解決風洞的邊界問題。在遷移完成后，域的邊角處某些 fi 會變得未知。原理圖（如下圖所示）顯示哪些 fi 對于其各自的邊和角未知。未知數(shù)由虛線箭頭表示：

為了理解未知的 fi 如何計算，讓我們考慮一下風洞的頂壁。對于此邊，f6, f7, f8是未知的。我們令 ：

是二維未知向量，而

 

 和 fieq 是平衡分布函數(shù)，定義為：     

 

Ho, Chang, Lin 和 Lin在他們的論文中詳細介紹了這種方法。此運算（實際上是一個高度可并行的運算）可以在 Mathematica 中有效地完成 ：

請注意，fieq 完全由速度和密度決定。將其代入速度和密度方程式，我們得到：

...其中 UBC, VBC 是用戶為邊界指定的速度。生成的這些方程是線性的。有三個未知數(shù) {ρ, Qx, Qy}，并且有三個等式，可以很容易地使用 LinearSolve 或 Solve 求解。同樣的步驟適用于所有邊和角。由于這是一個小的3 × 3系統(tǒng)，我們可以符號式地預先計算缺失分布的表達式。因此，對于頂壁，我們將有：

一旦 {ρ, Qx, Qy} 已知，將這些值代入即可獲取未知數(shù)：

在處理溢出邊界條件時，我們實際上是試圖跨邊界施加 0-梯度條件，即

 

。最簡單的一件事就是使用后面一個網(wǎng)格的速度，并用它來計算邊界上的fi。然而這將導致不準確的結(jié)果，并且在很多情況下，這些結(jié)果也將不穩(wěn)定。

 

 

考慮右壁，如下圖所示：

遷移步驟之后，分布 f4, f5, f6 未知。為了對分布函數(shù) fi 施加溢出條件，我們使用以下關(guān)系：

 

…… 其中

 uk(t) 

是上一個時間步在網(wǎng)格點(j, k), j = 1, 2, … 的速度，M 和 fi,N(t) 是上一個時間步長的分布。標量項u*必須為正數(shù)。一個好的候選是邊界處的正常速度，所以對于左右壁，有

 

，對于上下壁，

 

 

。如果u*是0，它可以被特征格子速度所取代。請注意，未知分布會根據(jù)哪面壁具有溢出條件而變化。例如，如果它是頂壁，則f6, f7, f8未知。Yang的文章（https://www.sciencedirect.com/science/article/pii/S0898122112006736）詳細介紹了這種溢出處理方法。 

 

 

而第二種方法只需做以下操作：

此方法比第一種方法簡單得多，并且基于對每個分布函數(shù)應(yīng)用二階逆向差分公式。

 

一旦遷移完成并施加邊界條件后，將進行"碰撞"（回憶一下基于規(guī)則的方法）。此步驟基本上是找出整體平均粒子將如何相互交互。此步驟有點復雜，但使用 BKG 近似，碰撞步驟可以寫作：

……其中是遷移步驟和邊界調(diào)整步驟后計算的密度和速度。vLBM 是格子玻爾茲曼域中的運動粘度，τ 是松弛參數(shù)。這個項非常重要，這里詳細介紹一下。

 

這在 Mathematica 中是兩行非常簡單的代碼：

從視覺角度來看，碰撞后分布函數(shù)根據(jù)前面的公式進行調(diào)整，我們得到了新的分布：

請注意箭頭長度的變化。就是這個！這就是運行格子玻爾茲曼模擬所需要的一切。

 

使用雷諾數(shù)將模擬變?yōu)楝F(xiàn)實

 

那么，如何模擬傳說中的納維-斯托克斯方程？為了回答這個問題，我們必須進行大量的數(shù)學和多尺度分析，但簡而言之，feq的形式?jīng)Q定了LBM模擬的宏觀方程。因此，使用適當?shù)钠胶夂瘮?shù)可以模擬一大堆偏微分方程。讀者有興趣的話，可以去網(wǎng)上找到很多人們使用LBM方法進行模擬的精彩范例。

 

我們已經(jīng)了解了LBM的基本機制，下一個顯而易見的問題是：在基于格子的系統(tǒng)中執(zhí)行的模擬如何轉(zhuǎn)換到物理世界中來？這是通過匹配格子系統(tǒng)和物理系統(tǒng)的非維參數(shù)來完成的. 在我們的例子中， 非維度參數(shù)只有一個 ：雷諾數(shù)。雷諾數(shù) （Rey） 的定義如下：

 

…其中Uphy是一個特征速度，Lphy是特征長度，vphy是物理域中的運動粘度。為了模擬流體，用戶需要指定雷諾數(shù)、特征速度和特征長度。這三條信息確定了運動粘度，隨后確定了松弛參數(shù)τ。使用這些信息和與 LBM 關(guān)聯(lián)的基礎(chǔ)方程，計算模擬的內(nèi)部參數(shù)。

 

 

特征格子速度 ULBM 不能超過（此數(shù)字計算為格子域中的聲速）。格子速度必須明顯低于此值，才能正確模擬不壓縮性。通常，格子速度為ULBM= 0.1。同樣，特征格子長度LLBM表示格子域中用于表示物理域中特征長度的點數(shù)。LLBM將是一個整數(shù)，并且通常由用戶定義。

 

讓我們通過一個例子看一下如何將格子模擬與物理模擬聯(lián)系起來。讓我們假設(shè) Rey = 100，Uphy = 1，Lphy = 1，風洞的物理尺寸為Lw = 15Lphy，Hw = 5Lphy. 對于格子域，我們假設(shè) ULBM = 0.1，LLBM = 10。格子風洞尺寸為LW,LBM = 15 × 10 = 150，Hw,LBM = 50。格子域中的粘度為：

這意味著 BGK 模型中的松弛參數(shù)τ為：

現(xiàn)在，我們擁有運行模擬所需的所有數(shù)量。如果我們現(xiàn)在讓模擬中的每個晶格時間步長為δtLBM = 1，那么我們需要知道δtphy是什么。這是通過等同粘度來完成，并且通過以下條件給出：

因此，如果我們想要運行的模擬 t = 100（時間單位），那么在晶格域中，我們將迭代

 

步。

 

 

雷諾數(shù)是一個顯著的非維參數(shù)。雷諾數(shù)不是指定我們模擬的流體及其速度和維度，而是將它們聯(lián)系在一起。也就是說，對于長度、速度和流體介質(zhì)截然不同的兩個系統(tǒng)，只要雷諾數(shù)保持一致，這兩個流體將具有相同的性能。

 

向風洞中添加物體

 

現(xiàn)在我們來討論下如何將物體放到風洞中。一種方法是將原始物體離散化成鋸齒狀，并將其與網(wǎng)格對齊，然后在每個步長的邊和角上施加無滑動邊界條件：

 

這種方法并不理想，因為它扭曲了原始物體。如要更好的表示物體，需要使用極其精細的網(wǎng)格，但這讓計算成本昂貴并且造成不必要的浪費。另外，尖銳的邊角會使流體產(chǎn)生不希望的行為。第二種方法是將物體浸入到網(wǎng)格中。物體的邊界被離散化，并浸入至域中：

 

 

對具體物體的邊界進行離散，可以使用內(nèi)置函數(shù)BoundaryDiscretizeRegion 輕松實現(xiàn)。我們可以指定 Disk 或 Circle 來生成一組點表示圓形物體的離散化版本：

 

這種將物體離散并放入網(wǎng)格的方法稱為浸入邊界法（IBM）。一旦離散化的物體被浸入，就需要模擬該物體對流體的影響，以確定其遵循邊界速度條件。一種使流體"感受"到浸入物體存在的方法是稱為直接施力方法。這種方法通過對演化方程添加施力項Fi ，對格子BGK模型進行了修正：

……其中為物體邊界在網(wǎng)格點上產(chǎn)生的矯正力。計算速度的公式現(xiàn)在修正為：

 

修正力的計算公式為：

 

 

……其中是物體的邊界條件，是不存在物體的情況下在物體邊界處的速度，是增量函數(shù)的近似值， 是物體邊界點的位置，通常稱為拉格朗日邊界點。有幾種選擇可以使用。我們將使用以下緊湊支持的功能：

 

 

這種近似也稱為柔化函數(shù)內(nèi)核，可以使用分段函數(shù)進行定義：

 

 

 

 

由此得到二維函數(shù) δ(x – X, y – Y)：

 ...其中dx，dy是縮放參數(shù)。對于拉格朗日點(Xi, Yi)，指定了δ函數(shù)。下面是以（1/2，1/2）為中心并按1縮放的δ函數(shù)的外觀 ：

讓我們通過一個示例來說明這個沉浸式邊界的概念，以及如何構(gòu)造函數(shù)并將其用于逼近函數(shù)。假設(shè)一個圓浸入矩形域中：

 

 

每個拉格朗日邊界點（藍色）都會影響一定半徑內(nèi)的網(wǎng)格點（垂直線和水平線的各個交點），如下圖所示：

為了獲得受每個拉格朗日點影響的網(wǎng)格點，我們使用了Nearest函數(shù)：

離散點的函數(shù)δ(x – X(s), y – Y(s))實際上成了一個矩陣：

該矩陣現(xiàn)在可以用于計算拉格朗日點的值。例如，讓我們假設(shè)底層網(wǎng)格上具有由h(x, y) = Sin(x + y)定義的值，則拉格朗日點的值計算如下：

...其中D是δ的離散化，h(xj, yj)是在(xj, yj)處要計算的函數(shù)值：

我們可以將計算出的插值與實際值進行比較：

同樣，可以使用拉格朗日點上的函數(shù)值來計算網(wǎng)格上的函數(shù)值：

如您所見，由于δ函數(shù)具有緊支集，因此只有位于其影響半徑內(nèi)的網(wǎng)格點才會獲得插值。所有不在支撐半徑內(nèi)的網(wǎng)格點均為0。

 

因此，在這里提醒讀者，格子玻爾茲曼模擬中的一個步驟包括以下步驟：

執(zhí)行遷移步驟。

在邊界處調(diào)整分布函數(shù)。

執(zhí)行碰撞步驟。

計算速度。

計算物體在拉格朗日邊界點的速度。

對于物體的每個邊界點，計算在該點施加邊界條件所需的校正力。

使用在步驟6中獲得的力計算點陣網(wǎng)格點處的校正力。

執(zhí)行遷移和碰撞步驟，考慮到作用力。

計算密度和速度。

 

到此，我們對使用二維LBM運行風洞模擬所需的所有要素已經(jīng)介紹完畢。

 

示例

 

為了使此風洞易于使用，所有函數(shù)都放入了名為 WindTunnel2DLBM 的程序包(https://blog.wolfram.com/data/uploads/2019/10/WindTunnel2DLBM.zip) 中。它包含許多功能，并且用戶可以輕松對問題進行設(shè)置。建議有興趣的用戶仔細閱讀程序包文檔，以了解詳細信息。這里將重點介紹各種示例并展示我們的計算型風洞設(shè)置的靈活性。

 

第一個示例是風洞中的流體。這也許是最簡單的情況。域及其相關(guān)邊界條件的示意圖如下所示：

 

在這種情況下，該問題只有一個長度尺度：風洞的高度。因此，這成為我們的特征長度尺度。在這種情況下，特征速度是來自進口的最大速度，該速度設(shè)置為1。剩下的就是指定要執(zhí)行模擬的雷諾數(shù)。這也是由用戶定義。讓我們將風洞的長度設(shè)為6個單位，從（0,6），將高度設(shè)為2個單位，從（-1,1）。現(xiàn)在，我們設(shè)置模擬：

 

請注意，我們在此處未提供任何邊界條件信息。這是因為風洞默認氣流在通道情況下，因此自動施加所有邊界條件。我們只需要指定特征信息和風洞的尺寸即可。

 

使用固定的時間步長執(zhí)行模擬。時間步長在內(nèi)部計算，可以從以下屬性訪問：

現(xiàn)在讓我們以5個時間單位運行模擬：

我們可以在模擬的最后一步查詢數(shù)據(jù)：

解決方案可以通過多種方式可視化。對于2D模擬，流線圖可以顯示一些有用的信息，讓我們可視化流線圖：

請注意，流線并非完全平行（存在一些偏差）。要了解原因，讓我們看一下x各個位置處速度場的u分量的分布：

這表明速度具有空間依賴性。對于此特定問題，我們應(yīng)該期望流量達到穩(wěn)定狀態(tài)，即流量不應(yīng)隨時間變化。讓我們以另外20個時間單位運行模擬，并查看速度曲線：

 

我們看到在各個x位置的速度分布幾乎彼此相同。這表明流確實達到了穩(wěn)定狀態(tài)。

 

 

箱內(nèi)流體問題 

 

現(xiàn)在，讓我們看一下經(jīng)典的"箱內(nèi)流體"問題。這是域的示意圖和邊界條件信息：

頂壁以水平速度1（長度單位/時間單位）移動，而其他所有壁都是固定防滑壁。箱內(nèi)的圓圈表示可能發(fā)生的流體行為。當頂壁移動時，壁會拖動下方流體，從而導致流體旋轉(zhuǎn)-這是原理圖中的大圓圈，它代表了渦旋。如果主旋渦有足夠的強度，則可能引發(fā)較小的次級旋渦。我們假設(shè)渦旋強度與雷諾數(shù)有關(guān)。讓我們看看以雷諾數(shù)100運行模擬會發(fā)生什么：

現(xiàn)在，我們迭代50個時間單位：

可視化結(jié)果表明，在中間附近形成了一個主旋渦，而箱的右下方形成了一個較小的次級渦旋：

如果雷諾數(shù)增加，則次級渦流會變得更強，并且額外的渦流會在拐角開始形成。讓我們看一下雷諾數(shù)為1,000的情況：

我們再次迭代50個時間單位：

可視化結(jié)果：

注意，主旋渦已移近中心。從它的外觀來看，已經(jīng)足夠強大到在域的左下角和右下角形成次級渦旋。

 

現(xiàn)在讓我們看看如果在"高"箱而不是正方形箱上進行仿真會發(fā)生什么。邊界條件保持不變，但域沿y方向變化：

運行模擬并使用 ProgressIndicator 跟蹤進度。此模擬將花費幾分鐘：

可視化流線：

在高箱場景中，在頂壁附近形成了一個主旋渦，而該旋渦又在其下方創(chuàng)建了另一個旋渦。如果第二個渦旋強度足夠大，它將在箱的底角產(chǎn)生渦旋。

 

我們已經(jīng)看到風洞為我們提供的靈活性。現(xiàn)在讓我們在風洞中放置一個物體并觀察氣流的行為。對于此示例，使用圓形物體：

這與通道中的流體相同（我們的第一個示例），但物體是放置在通道中。注意，現(xiàn)在存在兩個長度尺度d和H。特征長度的選擇盡管是任意的，但必須與流體的物理特性相聯(lián)系。在此示例中，如果要增大或減小物體的大小，則可以預期其后面的流型將發(fā)生變化。因此，自然的選擇是使用d作為特征長度。

 

讓我們將圓柱體放置在域的（3,0）處。令圓柱體的大小為1個長度單位。也就是特征尺度將為1。設(shè)域大小為 (0, 15) × (–2, 2)。該物體用 ParametricRegion指定：

在開始模擬之前可視化隧道是一個好主意，以確保物體位于正確的位置：

讓我們模擬10個時間單位的流體：

這里有兩點要注意：圓柱體后面的對稱渦流對和圓柱體內(nèi)的流體。一個特寫鏡頭表明，圓柱內(nèi)也有一些流型：

這種行為的發(fā)生是因為我們正在使用IBM。如前所述，IBM計算一組要施加到網(wǎng)格點上的力，以使表示該表面的表面速度為0。它不指定圓柱體內(nèi)應(yīng)該發(fā)生什么。因此，作為不可壓縮的流體，在圓柱體內(nèi)也存在流型。重要的是，物體邊界處的速度為0（無滑動）。

 

現(xiàn)在讓我們繼續(xù)迭代30個時間單位，看看圓柱體后面的流型發(fā)生了什么。有時，查看另一個稱為渦量（vorticity）的變量可能會幫助我們更好地了解正在發(fā)生的事情：

設(shè)置等值線的配色方案：

可視化渦量：

現(xiàn)在，我們注意到對稱性已被破壞，并被"波浪"行為所取代。旋渦清楚地表明了波浪狀的行為。我們在這里看到的被稱為圓柱體帶來的不穩(wěn)定性。這種不穩(wěn)定性繼續(xù)擴大，最終在圓柱體后方開始形成渦流。這種現(xiàn)象稱為"渦旋脫落"。在圓柱體表面產(chǎn)生了一個剪切層，該剪切層被帶到下游。

 

該渦旋脫落還取決于雷諾數(shù)。如果數(shù)值足夠小，則不會有任何脫落。但是，如果雷諾數(shù)達到100 - 150左右，則會觀察到脫落。為了正確地觀察這種現(xiàn)象，最好能看到這種流動的時間演變。第一步，通過定特征項和義隧道中的物體來設(shè)置問題：

為了產(chǎn)生渦量的時間演變，我們將在每個時間單位提取解并生成一系列圖形：

運行模擬可以清楚地看到，在圓柱體的背面形成了兩個渦流，剪切層緩慢地受到擾動，然后振幅逐漸增大，最終分裂成渦流脫落：

觀察移動物體引起的擾動

 

在下一個示例中，我們將采用沉浸邊界處理法，將一個圓形儲罐"浸入"風洞內(nèi)部。儲罐邊界處的速度為0。在此儲罐內(nèi)，我們將浸入一個橢圓形物體。該物體放置在儲罐壁附近，并在儲罐邊界遵循圓形路徑。格子玻爾茲曼法對于浸沒邊界的的靈活性使我們對移動物體有極大的靈活性。目的是研究當該物體移動通過靜止流體時會產(chǎn)生什么樣的擾動。

 

通過定義特征項來設(shè)置問題。在這種情況下，將以雷諾數(shù)400進行仿真。特征長度和速度為1。隧道中有兩個物體。第一個是固定的大型圓形儲罐。第二個是將在此儲罐內(nèi)移動的橢圓形物體：

檢查潛在問題的幾何形狀永遠是一個好主意。我們可以通過在初始化時簡單地提取離散化對象來實現(xiàn)。我們可以看到一切各歸其位：

和剛才一樣，我們將研究流體的渦量等值線。讓我們首先定義配色方案和將要繪制的等值線的水平：

現(xiàn)在模擬以60個時間單位運行：

運行流體擾動的時間演變表明，在達到更均勻的圓形擾動之前，最初在儲罐內(nèi)形成了非常漂亮的幾何圖案：

在彎曲和有障礙物的管道中的流體

出于好奇（和好玩），在下面的幾何體中會有怎樣的流型？

流體從右端進入管道，然后向上移動，然后從左端排出。塞子在右端的作用是什么？它會如何影響排水？

 

使用我們當前的設(shè)置，問題非常容易弄清楚。同樣，我們只是將管道和其中的障礙物浸入風洞中。風洞的左邊界、右邊界和頂邊界被指定為零速度。底部邊界的流出條件為-1 <= x <=-0.7，0速度條件為-0.7 <= x <= 0.7，拋物線速度曲線為0.7 <= x <= 1：

運行40個時間單位：

讓我們繪制渦量：

我們看到障礙物/塞子會產(chǎn)生渦流脫落，并沿著管道向下移動。讓我們看一下y = 0處的速度：

如果比較出口（左側(cè)）和入口（右側(cè)）之間的速度曲線，會看到出口速度幾乎是入口速度的一半。這足以表明，塞子使管道左端排出的液體減少，這和我們的預期一樣。

 

模擬翼型上的流體

 

作為最后一個例子，讓我們看一下機翼周圍的流體。翼型或稱翼剖面，是使得飛機升空的基本要素。機翼的類型很多，但我們將重點介紹一種簡單的機翼，由參數(shù)方程

 

表示。參數(shù)a控制翼型的厚度，參數(shù)b控制翼型的曲率：

為了使飛機"升空"，即離開地面，機翼頂面的壓力分布應(yīng)低于底面的壓力分布。這種壓力差會導致機翼（以及與其相連的所有物體）向上抬升。通過讓風以相當高的速度吹過其表面來實現(xiàn)該壓力差。第二個考慮因素是機翼通常需要傾斜以具有"攻角"。這樣做可以確保更大的提升。我們還將給機翼一個-10\[Degree]的攻角。該模擬將在雷諾數(shù)為1,000的情況下運行?，F(xiàn)在要指出的是，雷諾數(shù)為1000是一個很小的值。小型飛機的典型雷諾數(shù)約為100萬。由于網(wǎng)格尺寸較大，無法在筆記本電腦上進行全尺度模擬。但是，即使雷諾數(shù)只有1,000，我們也能夠很好地理解其潛在的動力學原理。在這個例子里，均勻的流體填充從左側(cè)進入。頂部和底部隧道邊界設(shè)置為周期性，右側(cè)邊界設(shè)置為流出。此處的特征長度將是機翼的厚度：

在開始模擬之前，提取離散化的對象，并檢查其是否位于風洞內(nèi)的適當位置：

注意到有大量的網(wǎng)格點。這是因為我們允許20個晶格點來解決薄機翼的問題?，F(xiàn)在，我們以10個時間單位運行模擬，這需要一點時間才能完成，原因是：（a）分辨率（即運行此模擬所需的網(wǎng)格點數(shù)）非常大 (800×200)；（b）要完成模擬，必須執(zhí)行20,000次迭代：

開始迭代過程：

讓我們首先看一下渦量圖：

就像鈍體的情形一樣，我們看到了渦旋脫落。對于機翼而言，這實際上不是理想的屬性。理想情況下，我們希望流體能夠擁抱表面。當氣流分開時（如您在機翼的頂面所看到的），壓降無法正確實現(xiàn)，并且機翼將無法產(chǎn)生升力。

 

現(xiàn)在讓我們看看壓力。除了繪制壓力外，我們還將繪制一個無量綱的參數(shù)，稱為壓力系數(shù)，由 Cp = 2(p – p∞)/(ρLBM U2LBM) 定義，其中 p∞ 是遠上游的壓力。我們對物體表面的壓力感興趣：

您會注意到這里有兩條直線。下方的直線代表頂面上的壓力，而上方的線代表底面上的壓力。顯然，盡管與翼型有所分離，但我們?nèi)詴@得一些壓力差。我們還可以繪制壓力等值線，并在翼型附近對其可視化：

如果仔細查看配色方案，的確會發(fā)現(xiàn)頂面壓力小于底面壓力。因此，這種翼型也許有希望。我們剛剛探討的流體動力學特性被稱為伯努利原理，該原理已在航空（如我們在此處看到的）和汽車工程等領(lǐng)域得到應(yīng)用。

這不過是開始，還有更多示例等你去親自嘗試！我們這里提供另辟蹊徑，來用 Mathematica 來研究并模擬流體動力學問題。LBM與IBM結(jié)合使用，對研究和分析流體流動是一個很好的工具。借助 Mathematica 的內(nèi)置函數(shù)，實現(xiàn)數(shù)字風洞的組裝非常簡單。軟件包 WindTunnel2DLBM 幫助我探索了流體動力學領(lǐng)域中許多引人入勝的概念（以及令人嘆為觀止的可視化）。我希望您也能從中得到啟發(fā)，并深入研究流體流動現(xiàn)象。