摘要:本文提供多尺度氣液兩相層/湍流動(dòng)數(shù)值仿真解決方案。本文主要討論Level Set方法在各類多相流動(dòng)中的應(yīng)用情況，包括耦合兩相流與共軛傳熱，并且提供了基于VirtualFlow軟件進(jìn)行數(shù)值處理的各種工業(yè)算例。

?

1.簡介

多相流以各種形式出現(xiàn)在環(huán)境工程和工業(yè)工程中，其往往具有相間傳熱傳質(zhì)的特點(diǎn)，包括自然界或燃燒室中的蒸發(fā)分散相、蒸汽輪機(jī)中的冷凝、液體噴霧、碳?xì)浠衔镙斔?、可溶性氣體吸收等過程?？梢钥吹胶朔磻?yīng)堆的熱工水力學(xué)和油氣輸送可能是兩相流存在最多的兩個(gè)領(lǐng)域。計(jì)算熱工水力學(xué)已經(jīng)為發(fā)生在輕水反應(yīng)堆瞬態(tài)過程中的復(fù)雜流動(dòng)和熱/質(zhì)傳輸現(xiàn)象提供了越來越復(fù)雜且可靠的數(shù)值分析。本文旨在提供基于Level Set技術(shù)，針對(duì)多尺度多相層/湍流動(dòng)的數(shù)值方案。

?

2.多相流的性質(zhì)與形式

多相流在自然界以及工業(yè)界以各種形式出現(xiàn)，其分類及流動(dòng)性質(zhì)受系統(tǒng)中的流體性質(zhì)及其所占據(jù)的體積分?jǐn)?shù)控制。如果輸送相在載體相中被嚴(yán)重稀釋，即流體之間不進(jìn)行大量的動(dòng)量交換，則該多相流系統(tǒng)可以被視為“混合流”，輸運(yùn)相通過分子和湍流效應(yīng)擴(kuò)散到載體介質(zhì)中。在Boussinesq流中，不超過10%，例如雪崩、污染物運(yùn)動(dòng)或沉積物運(yùn)輸?shù)?。如果輸運(yùn)相的特征長度尺度足夠大，使得界面動(dòng)量交換機(jī)制變得重要，則該多相流系統(tǒng)可以被歸類為“互穿介質(zhì)”。這類問題包括非受限氣泡羽流和碳?xì)浠衔锕艿懒鲃?dòng)，適用于氣泡和攪拌流形式。此外，我們還可以定義“界面流”，各相被明顯且大尺度的界面分開。分層管流與水波都屬于界面流。而形狀或密度不同的小實(shí)體，如固體顆粒、液滴和微小的氣體氣泡等，通常形成一個(gè)“分散流”系統(tǒng)。

圖1所示的波浪破碎問題涉及不同的長度尺度，作用于不同的時(shí)間尺度。波浪在壓力和風(fēng)切應(yīng)力作用下，從總流中提取動(dòng)能而形成；波浪最終會(huì)破碎分解成小尺度流動(dòng)，并將動(dòng)能耗散為熱量。波峰翻卷后會(huì)形成白浪或微破碎層，作為一個(gè)混合兩相流區(qū)域，該區(qū)域的多相流體很難區(qū)分。湍流在流動(dòng)中無處不在，在剪切面的兩側(cè)都有湍流流動(dòng)。大尺度湍流在頂部由氣流與波浪夾帶的作用形成。核心液體的流動(dòng)也屬于湍流，具有寬帶湍流譜。小尺度湍流可能在波峰處演化。由破碎波峰俯沖運(yùn)動(dòng)產(chǎn)生的海沫液滴在空氣側(cè)湍流的作用下分散并沉積。海底沉積物的輸運(yùn)受到波浪破潰輸送湍流的影響。在相間傳質(zhì)方面，海水霧滴通過蒸發(fā)導(dǎo)致氣-水的通量不守恒，而微小的水波破碎過程則增強(qiáng)了CO2等可溶性氣體的吸收。

3.實(shí)際兩相流問題的現(xiàn)代計(jì)算多相流體力學(xué)（CMFD）技術(shù)

3.1從平均雙流體建模到界面跟蹤方法

多相流和傳熱的建模是一門復(fù)雜的學(xué)科，其通?；谝恍┙?jīng)驗(yàn)公式建立。分析建模已經(jīng)從集總參數(shù)和一維模型(例如：（平衡）均相流模型、相滑移模型、漂移通量模型和雙流體模型等)發(fā)展到準(zhǔn)多維子通道模型，再到三維計(jì)算多相流動(dòng)力學(xué)(CMFD)模型，以及與DNS和LES湍流模型等耦合。CMFD的核心是準(zhǔn)確識(shí)別界面動(dòng)力學(xué)，通過該界面動(dòng)力學(xué)可以定義流動(dòng)形式，并可以量化相關(guān)的相間質(zhì)/熱傳遞機(jī)制。

自由表面和界面流動(dòng)是指兩相流動(dòng)問題中涉及兩種或兩種以上的不相混流體被隨時(shí)間演變的尖銳界面分開。通常，當(dāng)界面一側(cè)的流體是在界面上施加切應(yīng)力的氣體時(shí)，后者被稱為自由表面。界面跟蹤/捕獲方式能夠定位自由表面，其不是通過拉格朗日方法來跟蹤界面(例如通過跟蹤位于界面上的標(biāo)記點(diǎn)等)，而是基于歐拉意義上(固定網(wǎng)格)，保持跟蹤場的演變來捕獲自由界面，例如Level Set函數(shù)，體積分?jǐn)?shù)場等。

目前學(xué)界主要有兩種計(jì)算方法來求解捕捉界面流動(dòng)，可以分為前端跟蹤/捕獲和接口跟蹤/捕獲。前端跟蹤是指顯式參數(shù)化界面并及時(shí)精確跟蹤表面上點(diǎn)的方法，包括Stokes或勢(shì)流的邊界積分方法。其中，在二維或軸對(duì)稱幾何的情況下，界面只是一條曲線，通常通過弧長和時(shí)間來參數(shù)化。在前端跟蹤中，隨著界面的演變，標(biāo)記在拉格朗日意義上被跟蹤，盡管底層流場可以在固定的歐拉網(wǎng)格上求解。前端跟蹤方法不是本文欲闡述的重點(diǎn)。

接口跟蹤方法作為第二種方法，是本文闡述并采用其對(duì)多相流動(dòng)進(jìn)行求解的重點(diǎn)。接口跟蹤方法通過隱式表示捕獲接口來處理此類問題，包括Level Set方法（該方法將界面視為在所有空間上定義的函數(shù)的水平表面），以及VOF方法，該方法通過跟蹤網(wǎng)格中每個(gè)計(jì)算單元相對(duì)于一個(gè)流體相的體積分?jǐn)?shù)來捕獲界面的位置。體積分?jǐn)?shù)為0或1的單元格不包含界面，而具有分?jǐn)?shù)值的單元格則包含界面。

?

3.2界面跟蹤方法

界面跟蹤方法用于預(yù)測(cè)需要精確界面識(shí)別的兩相流，例如氣泡，液滴或液體射流的破裂等。該方法的關(guān)鍵是使用了具有可變流體性質(zhì)以及描述表面張力的單流體守恒方程。接下來主要描述VirtualFlow中采用的Level Set方法。在不可壓縮流動(dòng)條件下，用單流體形式表示的耦合流體運(yùn)動(dòng)和傳熱方程采用如下形式：

其中ρ是密度，p是壓力，μ是粘度。方程(2)中的右端項(xiàng)表示表面張力及其接觸線壁面貢獻(xiàn)(即項(xiàng)c)。另外，n表示界面的法向量，κ表示表面曲率，γ表示流體的表面張力系數(shù)，δ表示以界面為中心的光滑狄拉克函數(shù)。方程(3)中T為溫度，Cp為流體的熱容，λ為導(dǎo)熱系數(shù)，Q'''為體積熱源。

在VirtualFlow采用的Level Set方法中，互不摻混的各相之間的界面用連續(xù)函數(shù)φ表示，φ表示當(dāng)?shù)氐浇缑娴木嚯x，在界面上設(shè)為零，在一側(cè)為正，在另一側(cè)為負(fù)。通過這種方法，可以識(shí)別兩種流體，從而將物理界面的位置與連續(xù)函數(shù)相關(guān)聯(lián)。方程(2)中的流體性質(zhì)、體力和表面力都局部依賴于Level Set函數(shù)，其演化方程(存在相變時(shí))為:

式中，為傳質(zhì)速率，其既可以通過界面上的能量跳躍直接確定，也可以使用傳熱相關(guān)性進(jìn)行建模。如果傳質(zhì)率被迫僅在三線處起作用，則方程(2)中的源項(xiàng)僅應(yīng)用于包含三線的單元格。如果單位長度的傳質(zhì)率已知，則根據(jù)接觸線長度和單元面積將其轉(zhuǎn)換為單位面積的傳質(zhì)率。而流體性質(zhì)，如密度，粘度，熱容量和導(dǎo)熱系數(shù)等則基于Level Set函數(shù)進(jìn)行局部更新，并使用平滑的Heaviside函數(shù)在界面上進(jìn)行平滑:

在實(shí)際算例中，Level Set函數(shù)在按方程(4)經(jīng)過一次迭代步后，不再是距離界面的帶符號(hào)距離。為了恢復(fù)其在界面附近的正確分布，將re-distancing方程解為穩(wěn)態(tài)形式:

其中sgn(x)是Signum函數(shù)。采用非振蕩三階WENO格式，在(4)的每個(gè)迭代步后求解上式。采用三階龍格-庫塔顯式時(shí)間積分格式求解Navier-Stokes方程和Level Set方程。對(duì)流通量采用三階Quick格式進(jìn)行離散。采用二階中心格式差分?jǐn)U散通量。

3.3 VirtualFlow多相流軟件

VirtualFlow是一個(gè)基于求解多流體Navier-Stokes方程的多物理場、有限體積計(jì)算流體力學(xué)軟件。軟件使用結(jié)構(gòu)化網(wǎng)格，但允許將多個(gè)塊設(shè)置在一起。針對(duì)多塊問題，采用了基于MPI并行算法。網(wǎng)格排列是并置的，因此可以更容易地處理曲線傾斜網(wǎng)格。求解器是基于壓力的(投影型)，使用Karki-Patankar技術(shù)對(duì)低馬赫數(shù)可壓縮流進(jìn)行校正?？梢圆捎酶唠A時(shí)間推進(jìn)和對(duì)流格式;空間中的三階單調(diào)格式。多相流是用層流和湍流的界面跟蹤技術(shù)來解決的，在層流和湍流中，流動(dòng)被認(rèn)為是一種具有可變物質(zhì)特性的流體，這些特性根據(jù)被流動(dòng)平流時(shí)的顏色函數(shù)而變化，從液相中識(shí)別出氣體流動(dòng)區(qū)域。具體來說，Level-Set和VOF方法(ITM)都可以用于跟蹤不斷變化的界面。

?

3.4浸入式邊界/表面技術(shù)

浸沒邊界(IB)技術(shù)是在90年代后期發(fā)展起來的，用于模擬與實(shí)體邊界相互作用的流動(dòng)。浸沒表面技術(shù)(IST)的靈感來自于兩相流(VOF和Level Set)的界面跟蹤技術(shù)，其中自由表面由一個(gè)雙曲對(duì)流方程來描述。在IST中，固體被描述為第二相，具有自己的熱機(jī)械性能。該技術(shù)與Peskin的浸入邊界法有很大的不同，其主要優(yōu)點(diǎn)是可以更好地解決共軛傳熱問題。如圖2所示，流體繞流過一圓形物體，圓形實(shí)體浸入到由笛卡爾網(wǎng)格覆蓋的立方體網(wǎng)格中。實(shí)體由其外部邊界定義，使用實(shí)體水平集函數(shù)φs表示。與對(duì)多相流的分界面的處理一樣，φs是到表面的距離；在分界面上是零，在流體中為負(fù)，在固體中為正。固體表面粘滯剪切的處理方法與所有CFD程序的處理方法基本相同，其中壁面指數(shù)以邊界條件的形式表達(dá)。

?

4.計(jì)算算例

4.1 電子PCB板的冷卻

工程上往往對(duì)風(fēng)扇對(duì)電子電路板PCB的冷卻效果感興趣，特別是對(duì)不同特性的風(fēng)扇放置在不同的位置對(duì)PCB的冷卻效果的影響。CMFD有助于提供在特定設(shè)置下可能發(fā)生的情況的數(shù)值計(jì)算結(jié)果。VirtualFlow在該領(lǐng)域旨在幫助重新設(shè)計(jì)工程操作。得益于浸入式曲面模塊，VirtualFlow生成復(fù)雜三維算例設(shè)置(網(wǎng)格、區(qū)域、風(fēng)扇尺寸、材料屬性、邊界條件)僅約需要一個(gè)小時(shí)。

圖3反映了EST模塊在簡單笛卡爾網(wǎng)格中映射復(fù)雜形狀組件的能力，并在必要的地方(圍繞組件)進(jìn)行了細(xì)化。湍流的k-ε模型結(jié)合壁面函數(shù)進(jìn)行穩(wěn)態(tài)解析。網(wǎng)格由125000個(gè)單元格組成。如第3.3節(jié)所述，只采用高階方案。該圖描述了風(fēng)扇產(chǎn)生的冷卻流對(duì)從組件到堆芯流的傳熱的影響。在這個(gè)模塊中，可以自動(dòng)設(shè)置風(fēng)扇的特性，例如指定風(fēng)扇的數(shù)量、輪轂直徑、轉(zhuǎn)速和每個(gè)風(fēng)扇的平均流量。薄板元素被人為加厚，以使它們能夠由網(wǎng)格表示(至少3個(gè)單元應(yīng)該覆蓋內(nèi)部固體區(qū)域)。每個(gè)部件都可以指定其熱機(jī)械性能和內(nèi)部體積熱源。

4.2換熱器中的兩相流換熱

這里討論了耦合兩相流體流動(dòng)和共軛傳熱問題，測(cè)試案例是處理熱管理在熱交換器中的應(yīng)用，特別是在室內(nèi)空調(diào)大型熱泵中的使用。本驗(yàn)證算例將一根或一組有限壁厚的管道放置在自由湍流中，代表空氣冷卻劑。在管道內(nèi)部，熱水通過重力或外部泵送的作用流動(dòng)。算例比較了管道完全填充或部分填充的兩種極端情況。傳熱將由冷卻劑空氣向堆芯流動(dòng)觸發(fā)，通過三循環(huán)機(jī)制：在水中對(duì)流，通過管道固體壁傳導(dǎo)，最后從管道表面向堆芯外對(duì)流流動(dòng)。

采用IST方法對(duì)裝置(單管和組管)進(jìn)行網(wǎng)格劃分，如圖4所示。單個(gè)管道使用了包含98x20x66個(gè)單元格的笛卡爾網(wǎng)格，五個(gè)管道組使用了包含138x20x130個(gè)單元格的笛卡爾網(wǎng)格。用5個(gè)細(xì)胞覆蓋壁管。如上所述，所有方程均采用Quick格式。Level Set用于跟蹤管道內(nèi)的水面(當(dāng)它被部分填充時(shí);其中流入孔隙率設(shè)置為50%)。在穩(wěn)態(tài)條件下，使用k-ε湍流模型結(jié)合壁面函數(shù)僅對(duì)冷卻劑空氣進(jìn)行求解(假設(shè)入流湍流強(qiáng)度為5%)在管道內(nèi)流動(dòng)的水呈層流狀?；诠艿劳獍霃剑ㄌ卣鏖L度）的空氣流動(dòng)雷諾數(shù)Re=7000。

部分填充單管的模擬結(jié)果如圖5所示，其顯示了流入的冷卻劑，管道的自由表面，以及由冷卻劑空氣熱對(duì)流而排出的熱輪廓。

圖6顯示了管道的二維截面(位于中心)的流動(dòng)情況，顯示了完全填充和部分填充管道中的溫度等高線。圖6上圖顯示，由冷卻劑外部強(qiáng)制的管道傳熱在管道周圍是對(duì)稱的，與預(yù)期結(jié)果相符。事實(shí)上，如果流動(dòng)雷諾數(shù)達(dá)到一定程度，會(huì)引發(fā)渦脫落，此時(shí)非定常模擬就會(huì)顯示出一些對(duì)稱性的喪失。傳熱的非對(duì)稱性描繪在圖6下圖中，其中熱量主要向底部被液體占據(jù)的區(qū)域傳遞。

?

4.3輸油管道中烴類段塞的形成

數(shù)值模擬算例建立完整三維計(jì)算域，如圖7所示。管道長6.3m，管徑0.14m。如左圖所示，采用多塊網(wǎng)格策略，用相鄰子域覆蓋域。這里使用傳統(tǒng)的BFC網(wǎng)格和IST網(wǎng)格。網(wǎng)格子域分布在12個(gè)CPU上，在Linux PC集群上并行執(zhí)行MPI。計(jì)算網(wǎng)格數(shù)量大約30萬個(gè)。

?該仿真基于用于界面跟蹤的Level Set方法和用于湍流模擬的LES方法的結(jié)合。MILES方法用于模擬亞網(wǎng)格尺度(SGS)擴(kuò)散效應(yīng)。入口流速為14m/s，含氣率為50%。圖8顯示兩個(gè)交叉流位置在一個(gè)時(shí)間步長的自由液面位移預(yù)測(cè)，并給出了速度分布情況。

段塞的形成是由Kelvin-Helmholtz (KH)不穩(wěn)定性的第一模態(tài)增長觸發(fā)的。在第一模態(tài)的尾跡中形成的第二模態(tài)，后與第一模態(tài)合并最終導(dǎo)致段塞密封。段塞流在距離入口下游6個(gè)管徑的下游形成。隨著時(shí)間推移，由液體形成的堵塞會(huì)導(dǎo)致更多的液體積聚從而進(jìn)一步隔離了氣體段塞。該機(jī)制不斷重復(fù)發(fā)生，導(dǎo)致沿管道形成連續(xù)段塞。圖9顯示了正負(fù)渦度水平，表明在段塞破裂前，湍流在段塞的下游已特別活躍。

4.4燃料束子通道分析

燃料束子通道流動(dòng)問題屬于非常具有挑戰(zhàn)性的熱工水力學(xué)問題。參考Sadatomi等數(shù)據(jù)，模擬垂直2·3桿束子通道水力平衡兩相流的流動(dòng)特性。實(shí)驗(yàn)以常壓常溫下的水和空氣為工質(zhì)，2·3棒束通道為試驗(yàn)通道，包含6根矩形陣列棒束和2種6個(gè)子通道，模擬沸水堆燃料棒束(圖10)。在水力平衡流動(dòng)條件下，獲得了各種單相和兩相流動(dòng)下沿各子通道軸線的流量分布數(shù)據(jù)和壓降。裝置尺寸如下:桿徑:d = 16mm;桿距:p = 20mm;節(jié)徑比:p/d = 1.25;間隙:S11 = S12 = S22 = 4mm;液壓直徑:Dh =14.3 mm;流過面積:A = 194 mm2。計(jì)算域的長度減小到L = 120 mm，管的長度減小到Lt = 20 mm。域的寬度設(shè)置為W = 30 mm。通道底部和頂部分別設(shè)置流入和流出邊界條件。網(wǎng)格由65x65x130個(gè)單元格組成。數(shù)值模擬在一個(gè)8CPU的Linux PC集群上并行運(yùn)行，耗時(shí)24小時(shí)。

在預(yù)測(cè)沸水堆燃料棒束冷卻劑的熱水力特性時(shí)，必須準(zhǔn)確地評(píng)估子通道之間的流體傳遞。兩相系統(tǒng)中的流體傳遞由三個(gè)獨(dú)立的部分組成：空隙漂移、導(dǎo)流橫流和湍流混合。與前文中的測(cè)試算例相同，這里的流動(dòng)也是多尺度、多特征的，包括流體流動(dòng)、傳熱(共軛)、界面流動(dòng)，同時(shí)，流動(dòng)中的相變過程也至關(guān)重要。這里模擬了兩種情況:入水速度UL = 0.1 m/s的層流和入水速度UL = 1 m/s的弱湍流；流入孔隙率均為50%。圖11與圖12分別顯示了層流于湍流情況下，管與氣液界面之間的熱分布。

?

圖11與圖12的右側(cè)兩幅圖片為層流與湍流情況下不同的兩個(gè)時(shí)刻下的流場情況，圖片描述了管表面的溫度分布情況。由于中心管與氣相暴露在更高的表面交換面積下，這些管有更高的熱脫除率，這就解釋了為什么那里的溫度比鄰近的區(qū)域低(也因?yàn)橥ǖ劳鈬膶?duì)稱邊界條件)。在湍流情況，從溫度分布（顏色分布）顯示出流動(dòng)的高度變形，這導(dǎo)致管道的隨機(jī)熱傳遞過程。

?

5.結(jié)論

VirtualFlow對(duì)多尺度、多流體的層流和湍流氣液兩相流動(dòng)有很好的數(shù)值仿真效果，其采用多相流最先進(jìn)數(shù)值方法。VirtualFlow的亮點(diǎn)在于基于歐拉方法的界面跟蹤方法，特別是將Level Set技術(shù)用于模擬各種流動(dòng)。同時(shí)本文給出了基于Virtual Flow處理的各種工業(yè)實(shí)際算例：電子電路板的熱管理、管式換熱器內(nèi)三維耦合兩相流與共軛換熱、圓管內(nèi)段塞的形成以及子通道流量分析等算例。測(cè)試算例評(píng)估了應(yīng)用于兩相湍流的界面跟蹤方法的預(yù)測(cè)性能，實(shí)際算例的計(jì)算結(jié)果足以證明對(duì)復(fù)雜界面兩相流問題的數(shù)值仿真已是基于VirtualFlow的CMFD解決方案所能解決的。同時(shí)，在計(jì)算效率方面，本文中的計(jì)算算例均可在Linux PC集群上進(jìn)行數(shù)值計(jì)算，并能夠在24小時(shí)內(nèi)得到計(jì)算結(jié)果。