1、驗證OpenFOAM中interFOAM求解器 1. OpenFOAM簡介OpenFOAM是Open Field Operation and Manipulation的簡稱，顧名思義，它實現的是(張量)場的運算和操作，實質上是一個應用于計算連續介質力學的C+類庫。這一開源軟件起源于八十年代后期英國帝國理工大學的Gosman教授的團隊。當時的CFD代碼普遍采用FORTRAN語言編寫，為了尋求一種更為強大靈活通用的模擬平臺，在以后的若干年里，他們利用上了C+語言的高級特性，采用更為有效的方式重新編寫了很多代碼。其中最初的一些C+類出自Charlie Hill的博士論文，到1993年由Henry W

2、eller和Hrvoje Jasak以及后來的帝國理工博士生們持續不斷的開發。這一軟件起初叫做 FOAM，直到2004年由Henry創立的OpenCFD公司以開放源代碼的形式公之于眾，更名為OpenFOAM，目前最新版本是OpenFOAM-1.7，另外還有一個Jasak引導的克羅地亞版本，目前的版本是OpenFOAM-1.5-dev，兩個版本在底層類庫和上層應用的涵蓋范圍上有許多差異。OpenFOAM的優越性表現在以下幾個方面：1. 是最早利用C+語言編寫而成的科學軟件包之一（其它的主流CFD軟件公司已經發布或正在開發新一代C+代碼）；2. 利用C+的運算符重載功能使得頂層代碼在對偏微分方程的

3、描述上相對簡單了許多，且可讀性強，這使得OpenFOAM看上去就像一種非常自然的適用于模擬物理問題的編程語言；3. 是最早采用多面體單元網格的通用CFD軟件包，而這個功能得以實現是源于對模擬對象采用分層描述的自然結果；4. 是目前發布于開源許可下的最強大的通用CFD軟件包。OpenFOAM為CFD領域研究工作者們提供了一個強大的開源研究平臺，基于OpenFOAM展開的學習和研究目前正在歐美國家的研究院校和技術公司如火如荼的進行：Chalmer University of Technology從2006年開始開設了針對OpenFOAM的博士課程，他領導的渦輪機小組開發出用于模擬空泡和回旋流的代碼

4、，在瑞典還有另外幾所大學也開設了OpenFOAM相關課程，Niklas Nordin和他的研究小組開發出強大的模擬燃燒的求解器(dieselFoam)，瑞典的幾家公司開始將它們應用于工業；2008年丹麥的工業界和研究院校發起了提升國家開源CFD軟件的倡議，八家公司和兩所院校(Aalborg University和Technical University of Denmark)聯合開發丹麥的開源CFD軟件平臺，他們選用OpenFOAM開展了諸多項目；在德國，許多研究院校已經開始用OpenFOAM來工作，他們幾乎每半年就組織一次OpenFOAM的用戶交流研討會；挪威的Risavika燃氣能源公司開

5、展了數值模擬傳熱多相流研究二氧化碳的捕捉和封存，該研究機構已經發表了他們的研究成果并給出了含傳熱過程的可壓縮VOF兩相流求解器(compressibleInterFoam)的算例；另外，我們也可以從CFD-Online論壇上看到OpenFOAM板塊的帖子數僅次于ANSYS(Fluent and CFX)，其用戶交流活躍程度非同一般。當然，基于OpenFOAM開展工作也面臨著一定程度的挑戰和困難，相比較商業CFD軟件，OpenFOAM缺乏詳盡的用戶使用手冊，相關學習資料也很有限，它對研究工作者提出了更高的要求，需要研究者更加積極主動地去探索和持續不斷地深入學習。CFD-Online的OpenFO

6、AM板塊是很好學習交流平臺，OpenFOAMwiki是很好的資料庫，對審視OpenFOAM的所有相關內容很有幫助，Jasak的克羅地亞網站提供了很多演講報告和OpenFOAM相關的博士論文，是深入學習的好資料，官方發布的UserGuide和ProgrammerGuide，盡管寫得遠遠不夠詳盡，初學者會覺得對讀程序沒有多大幫助，但是只要對OpenFOAM自定義化的語言一定程度的熟悉之后，再回過頭來看，會對OpenFOAM的認識有一個整體上的提升。2. 基本兩相流求解器interFOAM簡介2.1簡介兩相流(two fluid system)問題interFOAM是OpenFOAM中用來解決two

7、-fluid systems一類問題的求解器，屬于the flow of two immiscible fluids的現象有：水波、潰壩、水中氣泡等，two-fluid system按照界面的拓撲結構可分為三類：分離(segregated)、混合流(mixed)、分散流(dispersed)。比如一個盛有部分水的封閉容器，當低頻小幅度晃動時，界面完好，此時為第一種情況；當晃動頻率振幅增加到一定時，界面變得不穩定，界面部分破碎，致使液體包進了氣泡，此屬第二類情況；當容器劇烈晃動時，大量空氣以懸浮的氣泡形式存在液體中，為第三類情況。interFOAM中沒有引入湍流模型，因此它主要針對第一二類問題，

8、所有的方法能夠模擬界面發生較大變形，破裂或融合等現象，但也不排除第三類情況。這類問題的關鍵技術在于要模擬出兩種流體的界面interface，由于在整個體系的運動過程中，界面往往是不斷變化的，整個流體域的運動相互耦合，從而需要在求解過程中不斷生成界面，要解決好這個問題，需要考慮以下幾點：1. 在離散網格中表示出界面2. 如何處理分布了兩種流體的計算單元3. 界面的運動4. 界面條件與運動方程的耦合總的來說，interFOAM中使用了VOF（volume of fluid）來處理這類問題。自由面或流體界面的模擬方法總的可歸為兩類：面法（假想界面上有粒子，高度函數，Level Set 法，貼面法等）

9、和體積法（假想流體中有粒子，體積分數法等）。interFOAM 使用了體積法(volume method)中的體積分數法（volume fraction）來表示計算域中兩種流體的分布，但是僅僅知道計算單元(cell)中的兩種流體體積比還是不夠的（圖1），還要能進一步表示出每個計算單元中的流體分界面，這需要良好的interface capturing method，以期望模擬出真實的界面。目前出現的技術有：線技術(line techniques)，施主受主法(donor-accepter method)，高階差分法(higher order differencing schemes)，inter

10、FOAM使用了基于施主受主法的VOF法。 圖1. 離散網格中的體積分數2.2 數學建模對于two-fluid system 這類問題，首先引入體積分數函數，定義表示一種流體的存在，表示另一種流體的存在，表示兩種流體的界面區域，因此（）表示了整個流場中兩種流體的分別，在建立運動方程時將兩種流體當作一種流體來處理，該流體的物理屬性表示為： （1） （2）OpenFOAM采用FVM法離散計算區域，函數在界面區域經歷非連續變化(a step function)，這導致成為分段連續函數，為了能夠使用連續介質力學方法建立數學模型并計算界面曲率，需要對函數進行光滑處理，定義流體的過渡區域厚度為，則有：（3）

11、滿足的輸運方程： （4）若用表示任何守恒物理量，它滿足控制體形式的守恒方程為：邊界處的對流項邊界處的擴散項源項邊界處的源項 （5）方程（5）的微分形式為： （6）取，表示單位體積質量，由于我們研究的系統中無化學反應或相變，因此無源無擴散，方程簡化為： （7）利用方程（1）（4）可進一步簡化為 （8）動量守恒時，取，表示單位體積動量，認為系統中無動量擴散，因此，源項由內部力和外部力組成，流體微元內部力相互抵消，微元表面存在應力，切向表現為粘性力，法向表現為壓力，把流體視為局部熱平衡的牛頓流，則應力張量表示為： （9）另外，液體在自由面上還受到表面張力的影響，是張力系數，表示自由面每增加單位面積所

12、需做的功，用表示界面曲率，表示界面的單位法向量，用表示表面張力沿界面法向的一個分量，此力的作用是平衡界面兩邊的壓力差，此力只在界面處存在，在非界面處其值為零。 （10）其中 （11） （12）因此，得到動量方程：由方程（8）（9）進一步可化為： （14）在interFOAM中，為了解決階段函數的遷移問題，Weller引入了額外的人工壓縮項： （15）其中是一個適于壓縮界面的速度場，這一項只對界面區域產生影響。函數在界面區域經歷非連續變化，這導致密度粘性成為分段連續函數，為了能夠使用連續介質力學方法建立數學模型并計算界面曲率，還需要采取對函數進行光滑處理等相關技術。2.3 數值方法OpenFOA

13、M采用基于任意非結構化網格的有限體積法，將流體區域分成有限數量的互不重疊的小控制體，然后將微分形式的控制方程在每個小控制體積分，保證方程局部和整體上的守恒性，然后用分段之和近似積分，并提供各種差分差值方案（Euler implicit or explicit, Crank-Nicolson, centred, upwind, TVD, NVD等）來估算各種微分運算（面法向梯度，法向梯度，拉普拉斯項，散度項，時間項，flux計算），來近似表現流體的遷移或擴散等行為，使用時要就具體問題加以選擇。在求解代數方程組時，也有多種求解控制器可供選擇（Linear Solver Control, PISO

14、and SIMPLE 等）。在用interface capturing方法處理界面的情況中，容易出現數值擴散，導致界面受污，為了克服此問題，常采用界面壓縮技術，在OpenFOAM中求解方程時，Weller引入了他自己的一種技術叫做interfaceCompression，來解決傳統的VOF界面壓縮法存在的基本問題。3潰壩實例驗證interFOAM利用interFOAM可以求解很大一部分屬于two-fluid systems 的問題，所有的方法上一章已經介紹過，各種具體的問題的差異性只在于體積分數，壓力，速度的初始條件和邊界條件的不同。3.1 潰壩實例一：(1). 初始條件和邊界條件圖2 初始時

15、刻的水柱潰壩的初始時刻如圖2所示，左右和底部是墻（無滑移），此處速度始終為零；因為模擬的是二維效應，所以模型的前后面設置成滑移條件（OpenFOAM 本是三維模擬）；頂部模擬的是開邊界，因此在此處設置空氣壓力為定值，速度法向梯度為零，兩種流體可通過此邊界離開計算區域，但只有空氣可經此流入計算域，因此速度指向外時，體積分數為零梯度，速度向內時，體積分數為定值；初始時刻，速度為零，邊界處壓力梯度為零，體積分數如圖 2 在水柱區域為 1，其他區域為零。為了便于后面與實驗結果比較，這里取得幾何模型與實驗中的一致。(2)在 OpenFOAM 中求解要用 OpenFOAM 模擬潰壩，需要模型的幾何信息，輸

16、運性質（密度粘性等），環境信息（重力場），這些信息分別保存在“constant”文件夾中；需要整個計算域初始時刻的體積分數，速度和壓力信息，它們分別保存在“0”文件夾中；需要控制方程各項的離散方案，離散后的代數方程的迭代方法和容許誤差等信息，以及求解時需要的時間設置和讀寫控制等信息，它們分別保存在“system”文件夾中。建立了上面所述的輸入信息之后，在 Linux 下的終端依次執行：“blockMesh”（網格生成）“setFields”（設置流場初始分布）“interFOAM”（求解）“paraFOAM”（后處理），最后可在 paraView 中得到可視化的結果，計算出的各個時刻的壓力速度

17、體積分數這些數據保存在以時間命名的文件夾里（正如初始時刻數據放在“0”文件夾中一樣）。針對 Weller的界面壓縮方法，需要對數值方案和壓縮系數進行多種嘗試和比較，從而找出最優方案。(3)結果分析潰壩實驗是用于驗證自由面流或兩相流數學模型合理性的經典實驗。這里將所得結果（圖 4）與經典實驗結果（圖 3）做定性的比較。圖 3 到圖 6 都是反映潰壩的時間演化過程的圖片。在此過程中，重力加速度作用于左邊的水柱向下向右運動以尋求可能的最低勢能。初始階段由慣性力支配，隨后粘性效應迅速增大。下面捕捉幾個主要現象加以比較和分析：1)t=0.2s：底部約 75%的部分被水覆蓋，與實驗基本一致；2)t=0.4

18、s:：運動在最前端的水流爬上對面的墻上，底部的水平自由面如圖3中一樣略微傾斜3)t=0.6s：底部的水面已經平行于底部，爬上對墻的水流因為重力開始下滑4)t=0.8s：下滑后的水流形成一股反向運動的波并包進了部分空氣，但是與圖 4 相比，可以觀察到空氣泡體積偏大，另外返回的波原本在水面上濺起的大量水珠，在數值模擬中未能捕捉到，以至原本分散到空氣中的這部分還保留在波頭中，使波頭看起來比實驗照片中要低。這說明模型還不能很好的模擬出更小尺度的破碎（break up）和擴散(dispersed)現象。5) t=1.0s：返回波爬升到左邊墻上，由于波頭在水面上運動的跳躍性，在爬上墻時包進了大量空氣，空氣

19、泡中和水面上的擴散（dispersed）現象同樣沒有捕捉到，而其他情況基本一致。圖 3 潰壩的實驗結果圖4 潰壩的數值模擬結果3.2 潰壩實例二：圖5. 有障礙物的潰壩初始幾何模型更有趣的潰壩現象是，在水柱不遠處增加了一個障礙物(如圖5)，相比上例此例只在幾何建模時略有差異，其他信息完全不變。模擬出的結果（如圖7）：1)t=0.1s：實驗中（圖 6）此刻底部流體前段已經觸到障礙物，而數值模擬結果顯示此刻還未到達障礙物；2)t=0.2s：波頭被障礙物擋住，一部分繞過障礙物從上方繼續前進；3)t=0.3s：波頭繼續朝對面墻運動，并已蓋住了障礙物的上邊緣；4)t=0.4s：繞過障礙物的波頭撞倒對面墻

20、上，將下方的空氣困在里面，情況跟實驗吻合；5) t=0.5s：此時障礙物上方的部分水流由于重力作用開始下墜，撞到墻上的水流也開始下滑，但下方困住的空氣的阻礙作用，可以看到墻上下墜的水流和流過障礙物上邊緣的水流，在它們下方將會形成一個氣泡（如圖 7 中黃色小框標記出的）；6)t=0.6s：被困住的大氣泡中的空氣已經沖破上面水層，還可以看到上一時刻中所標示的已經形成完全被水所包裹的氣泡，另外繞到障礙物后方的水流在撞到底部之后，又反濺起來（圖中已標示）。這些現象跟實驗吻合。同樣，整個過程中脫離水面擴散到空氣中的大量水珠未能被模擬出來。圖 6.潰壩撞擊障礙物的實驗結果（Koshizuka(1995)圖

21、7. 潰壩撞擊障礙物的數值結果3.3 小結：上面潰壩的數值模擬結果已經在一定程度上顯示出了interFOAM模擬這類two-fluid systems 問題的有效范圍，interFOAM對自由面變形破裂和融合的模擬效果還是不錯的，但對根據細微的水和空氣微粒相互擴散現象的模擬目前難以做到，需要更加細致的網格或模型。當然在OpenFOAM中，由于所用方法的不同，解決這類問題還可以選用其他求解器，比如rasInterFOAM，該求解器中引入了各種湍流模型和DPE模型(Dispersed Phase Element)，應該可以模擬現象更加復雜的問題。4.造波驗證interFOAM求解器4.1 使用速度

22、入口式造波(1) 數值造波理論由線性造波理論，對于平衡位置在原點、沖程為、角頻率為的活塞式造波機，其推板做簡諧運動的速度為： （16）在水深為的波浪水槽中距造波板x處的波面為 （17）式中滿足方程 （18）而為下面方程的第n個根 （19）式（17）中項為造波板產生的衰減立波，在離開造波板一定距離后很快就衰減掉，而第一項則是造波板所產生的波數為，頻率為的行進波。令式（17）中的x=0，則得到造波前的波面為 （20）其中，若令為普通造波機產生所需要的余弦波面的推板速度，則有 (21)(2) 在OpenFOAM中如何實現？靜網格法，通過設置入口處的速度、壓力、gamma值，達到造波的目的，通過在wa

23、veMakerVelocityFvPatchVectorField.C文件中加入# include “waveMakerVelocity.H”，來調用waveMakerVelocity.H文件。waveMakerVelocity.H文件內定義了多種類型的波l 相當于定義了新的速度邊界條件waveMakerVelocity，在case的定義中直接調用自定義的邊界條件l 有OpenFOAM自帶的邊界條件進行改寫，這里根據OpenFOAM-1.6/src/finiteVolume/fields/fvPatchFields/inletOutlet進行改寫。(3) 線性波的數值模擬結果建立一個長45米，高1.2米，水深為0.6米的二維水池模型，用速度入口式生產一個波高為0.02米，周期為1s的線性波。水池中生成的波面如圖8、圖9。生成的波面圖圖8 線性波的波面圖圖9 距離造波邊界處1個、2個、4個、6個波長處的波面時歷曲線4.2 如何在OpenFOAM求解器中加海綿層，實現消波的目的為了更有效地解決數值水槽末端的波浪反射問題，在VOF模型中采用了海綿層阻尼消波的方法。其原理是將某種形式的衰減系數加入到控制方程之中，使得在海綿層內部的流體滿