導(dǎo)語(yǔ)

 

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

 

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

創(chuàng)建我們的計(jì)算型風(fēng)洞

 

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

 

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

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

 

當(dāng)人們思考計(jì)算流體力學(xué)時(shí)，總會(huì)立刻想到著名的納維-斯托克斯方程。這組方程是決定流體流動(dòng)行為的控制方程。在這一點(diǎn)上，我們似乎要使用納維-斯托克斯方程來(lái)幫助我們建造一個(gè)風(fēng)洞。但事實(shí)證明，還有其他方法可以研究流體流動(dòng)的行為，而無(wú)需求解納維-斯托克斯方程。格子玻爾茲曼方法 （LBM）就是這些方法中的一個(gè)。

 

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

 

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

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

 

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

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

            
密度和速度的計(jì)算可以如下所示：

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

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

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

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

 

由于δtLBM = 1 且

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

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

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

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

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

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

是二維未知向量，而

 

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

 

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

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

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

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

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

 

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

 

 

考慮右壁，如下圖所示：

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

 

…… 其中

 uk(t) 

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

 

，對(duì)于上下壁，

 

 

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

 

 

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

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

 

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

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

 

這在 Mathematica 中是兩行非常簡(jiǎn)單的代碼：

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

請(qǐng)注意箭頭長(zhǎng)度的變化。就是這個(gè)！這就是運(yùn)行格子玻爾茲曼模擬所需要的一切。

 

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

 

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

 

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

 

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

 

 

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

 

讓我們通過(guò)一個(gè)例子看一下如何將格子模擬與物理模擬聯(lián)系起來(lái)。讓我們假設(shè) Rey = 100，Uphy = 1，Lphy = 1，風(fēng)洞的物理尺寸為L(zhǎng)w = 15Lphy，Hw = 5Lphy. 對(duì)于格子域，我們假設(shè) ULBM = 0.1，LLBM = 10。格子風(fēng)洞尺寸為L(zhǎng)W,LBM = 15 × 10 = 150，Hw,LBM = 50。格子域中的粘度為：

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

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

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

 

步。

 

 

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

 

向風(fēng)洞中添加物體

 

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

 

這種方法并不理想，因?yàn)樗で嗽嘉矬w。如要更好的表示物體，需要使用極其精細(xì)的網(wǎng)格，但這讓計(jì)算成本昂貴并且造成不必要的浪費(fèi)。另外，尖銳的邊角會(huì)使流體產(chǎn)生不希望的行為。第二種方法是將物體浸入到網(wǎng)格中。物體的邊界被離散化，并浸入至域中：

 

 

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

 

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

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

 

修正力的計(jì)算公式為：

 

 

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

 

 

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

 

 

 

 

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

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

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

 

 

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

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

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

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

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

我們可以將計(jì)算出的插值與實(shí)際值進(jìn)行比較：

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

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

 

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

執(zhí)行遷移步驟。

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

執(zhí)行碰撞步驟。

計(jì)算速度。

計(jì)算物體在拉格朗日邊界點(diǎn)的速度。

對(duì)于物體的每個(gè)邊界點(diǎn)，計(jì)算在該點(diǎn)施加邊界條件所需的校正力。

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

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

計(jì)算密度和速度。

 

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

 

示例

 

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

 

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

 

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

 

請(qǐng)注意，我們?cè)诖颂幬刺峁┤魏芜吔鐥l件信息。這是因?yàn)轱L(fēng)洞默認(rèn)氣流在通道情況下，因此自動(dòng)施加所有邊界條件。我們只需要指定特征信息和風(fēng)洞的尺寸即可。

 

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

現(xiàn)在讓我們以5個(gè)時(shí)間單位運(yùn)行模擬：

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

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

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

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

 

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

 

 

箱內(nèi)流體問(wèn)題 

 

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

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

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

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

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

我們?cè)俅蔚?0個(gè)時(shí)間單位：

可視化結(jié)果：

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

 

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

運(yùn)行模擬并使用 ProgressIndicator 跟蹤進(jìn)度。此模擬將花費(fèi)幾分鐘：

可視化流線：

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

 

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

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

 

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

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

讓我們模擬10個(gè)時(shí)間單位的流體：

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

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

 

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

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

可視化渦量：

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

 

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

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

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

觀察移動(dòng)物體引起的擾動(dòng)

 

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

 

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

檢查潛在問(wèn)題的幾何形狀永遠(yuǎn)是一個(gè)好主意。我們可以通過(guò)在初始化時(shí)簡(jiǎn)單地提取離散化對(duì)象來(lái)實(shí)現(xiàn)。我們可以看到一切各歸其位：

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

現(xiàn)在模擬以60個(gè)時(shí)間單位運(yùn)行：

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

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

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

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

 

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

運(yùn)行40個(gè)時(shí)間單位：

讓我們繪制渦量：

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

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

 

模擬翼型上的流體

 

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

 

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

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

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

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

開(kāi)始迭代過(guò)程：

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

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

 

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

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

如果仔細(xì)查看配色方案，的確會(huì)發(fā)現(xiàn)頂面壓力小于底面壓力。因此，這種翼型也許有希望。我們剛剛探討的流體動(dòng)力學(xué)特性被稱為伯努利原理，該原理已在航空（如我們?cè)诖颂幙吹降模┖推嚬こ痰阮I(lǐng)域得到應(yīng)用。

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