流體-固體的節(jié)點(diǎn)化兩相流建模方法
【技術(shù)領(lǐng)域】
[0001 ]本發(fā)明涉及流體-固體的節(jié)點(diǎn)化兩相流建模方法。
【背景技術(shù)】
[0002] 目前，重力驅(qū)動簡單兩相流固耦合傳熱問題在工程上主要有兩種求解方法。
[0003] 第一種方法用基于擴(kuò)散分界面的"相場"方法來求解兩相流區(qū)域的流場，并與固體 實(shí)時(shí)耦合。主要應(yīng)用于商業(yè)CFD(計(jì)算流體力學(xué)）軟件，如FLUENT。"相場"方法用序參數(shù)1 (x，t)、<i>b(x，t)對兩相a、b的空間分布進(jìn)行建模。在相a占據(jù)的區(qū)域<K(x，t)等于1，在相a不 存在的區(qū)域Φ a(X，t)等于0，在相a周圍的擴(kuò)散分界面內(nèi)Φ a(X，t)在1到0之間連續(xù)變化。在流 場中的任何位置，巾3(1，0與(^(1，〇之和都等于1。在擴(kuò)散分界面內(nèi)，流體的特征參數(shù)通過 對兩相各自的特征參數(shù)依據(jù)序參數(shù)進(jìn)行插值來獲得。所以，流場中流體的特征參數(shù)處處連 續(xù)，在整個(gè)求解區(qū)域內(nèi)各相的控制方程是統(tǒng)一的。這種方法適用范圍較廣，可以處理比較復(fù) 雜的相變兩相流問題，但計(jì)算效率較低，而且處理不互溶的簡單兩相流時(shí)有一定誤差。第二 種方法用節(jié)點(diǎn)化的方法對流體、固體進(jìn)行建模，求解由這些熱節(jié)點(diǎn)組成的節(jié)點(diǎn)熱網(wǎng)絡(luò)。主要 應(yīng)用于商業(yè)節(jié)點(diǎn)熱網(wǎng)絡(luò)軟件，如FIowmaster。在工程應(yīng)用中，這種算法將固體結(jié)構(gòu)、兩相流 體抽象成若干個(gè)具有一定溫度、質(zhì)量的節(jié)點(diǎn)?？紤]各節(jié)點(diǎn)間的對流、輻射、導(dǎo)熱作用，以及環(huán) 境輻射、油箱內(nèi)熱源，就可以對各節(jié)點(diǎn)列出熱平衡方程。對整個(gè)系統(tǒng)的熱平衡微分方程組進(jìn) 行瞬態(tài)迭代求解，就可以實(shí)現(xiàn)對整個(gè)流固耦合傳熱系統(tǒng)的仿真、求解。雖然這種方法具有高 的計(jì)算效率，但計(jì)算結(jié)果的質(zhì)量較差，無法給出準(zhǔn)確的溫度場。

【發(fā)明內(nèi)容】

[0004] 根據(jù)本發(fā)明的一個(gè)方面，提供了一種流體-固體的節(jié)點(diǎn)化兩相流建模方法，其特征 在于將兩相流體中的液相和氣相分別抽象為具有溫度、質(zhì)量、體積的節(jié)點(diǎn)，并包括：
[0005] A)計(jì)算當(dāng)前時(shí)刻兩相流體中的液相和氣相的體積，通過重力方向的體積積分計(jì)算 出兩相分界面；
[0006] B)使用對流換熱公式，計(jì)算流固耦合壁面的對流換熱熱流密度；
[0007] C)使用一個(gè)表面輻射換熱模型，計(jì)算流固耦合壁面的輻射熱流密度，并將輻射熱 流密度與步驟B)中求出的對流換熱熱流密度相加，獲得流固耦合壁面的總熱流密度；
[0008] D)計(jì)算下一時(shí)間步內(nèi)氣相和液相獲得的總對流換熱熱量，分別更新氣相和液相的 加權(quán)平均溫度；
[0009] E)根據(jù)流固耦合壁面的總熱流密度，計(jì)算流固耦合壁面的溫度梯度，進(jìn)行一個(gè)時(shí) 間步的固體瞬態(tài)傳熱計(jì)算；
[0010] F)下一時(shí)間步之后的時(shí)刻，重復(fù)A)至E)，直到到達(dá)瞬態(tài)傳熱終止時(shí)刻。
【附圖說明】
[0011]圖1是根據(jù)本發(fā)明的一個(gè)實(shí)施例的流體-固體的節(jié)點(diǎn)化兩相流建模方法的流程圖。
[0012] 圖2是一個(gè)具備隔熱層和內(nèi)置電子設(shè)備的典型飛機(jī)燃油箱模型的縱剖面示意圖， 用于檢驗(yàn)本發(fā)明的方法的可靠性。
[0013] 圖3A顯示了用根據(jù)本發(fā)明的流體-固體的節(jié)點(diǎn)化兩相流建模方法對圖2所示的飛 機(jī)燃油箱模型進(jìn)行的數(shù)值計(jì)算的結(jié)果；圖3B顯示了用Fluent軟件對圖2所示的飛機(jī)燃油箱 模型進(jìn)行的計(jì)算結(jié)果。
【具體實(shí)施方式】
[0014] 根據(jù)本發(fā)明的流體-固體的節(jié)點(diǎn)化兩相流建模方法，將兩相流體抽象成具有一定 溫度、質(zhì)量、體積的節(jié)點(diǎn)，用節(jié)點(diǎn)熱網(wǎng)絡(luò)的思想來處理，規(guī)避多相流CHH十算。同時(shí)，依然用數(shù) 值傳熱學(xué)的方法來求解固體區(qū)域的瞬態(tài)傳熱，保證固體區(qū)域溫度的計(jì)算結(jié)果質(zhì)量。由于計(jì) 算中將流體節(jié)點(diǎn)化處理，所以輻射模型選擇不依托流體區(qū)域的表面輻射換熱模型(S2S)。
[0015] 如圖1所示，根據(jù)本發(fā)明的一個(gè)實(shí)施例的流體-固體的節(jié)點(diǎn)化兩相流建模方法包 括：
[0016] 1)將液相、氣相分別抽象為具有溫度、質(zhì)量、體積的"節(jié)點(diǎn)"；
[0017] 2)計(jì)算當(dāng)前時(shí)刻兩相流體中氣相和液相的體積，通過重力方向的體積積分計(jì)算出 液面高度(兩相分界面）；（步驟101、102)
[0018] 3)使用對流換熱公式，計(jì)算流固耦合壁面對流換熱熱流密度；（步驟103)
[0019] 4)使用表面輻射換熱模型（S2S)模型，計(jì)算流固耦合壁面的輻射熱流密度（步驟 104)，并將輻射熱流密度與步驟103中求出的對流換熱熱流密度相加，獲得總熱流密度；（步 驟 105)
[0020] 5)計(jì)算下一時(shí)間步內(nèi)氣相和液相獲得的總對流換熱熱量，分別更新氣相和液相的 加權(quán)平均溫度；（步驟106、107)
[0021] 6)根據(jù)流固耦合壁面的總熱流密度計(jì)算流固耦合壁面溫度梯度，進(jìn)行一個(gè)時(shí)間步 的固體瞬態(tài)傳熱計(jì)算；（步驟108、109)
[0022] 7)在下一時(shí)刻（下一時(shí)間步之后的時(shí)刻），重復(fù)2)至6 )，直到到達(dá)瞬態(tài)傳熱終止時(shí) 刻。
[0023] 2.1兩相區(qū)域劃分
[0024]根據(jù)本發(fā)明的一個(gè)實(shí)施例，在如圖1中的步驟101、102所表示的處理中，為了確定 流固耦合壁面上的對流換熱系數(shù)，需要對兩相流體區(qū)域進(jìn)行劃分。
[0025] 假設(shè)液相質(zhì)量流速很小且液面的波動可以忽略，則兩相區(qū)域劃分可以用液相的液 面高度來代替。即首先根據(jù)進(jìn)/出口質(zhì)量流量獲得每一時(shí)刻液相體積，然后根據(jù)液相體積確 定液面高度，即兩相分界面。
[0026] 考慮到固體內(nèi)幾何模型往往比較復(fù)雜，液面高度可根據(jù)液相體積通過網(wǎng)格單元體 積分的方法更簡單、更準(zhǔn)確地獲得。具體方法包括：
[0027] 1)獲得所有體單元的中心點(diǎn)坐標(biāo)和體積，并根據(jù)中心點(diǎn)高度由低到高進(jìn)行排序。
[0028] 2)設(shè)定初始液面高度，將該高度以下所有體元的體積相加，即可近似得到初始液 相體積及質(zhì)量。
[0029] 3)根據(jù)進(jìn)/出口質(zhì)量流量獲得各時(shí)刻液相質(zhì)量，并結(jié)合液相平均溫度確定每一時(shí) 刻的液相體積。
[0030] 4)按中心點(diǎn)高度從低到高依次將對應(yīng)體單元的體積相加，直到第i體元時(shí)體積之 和首次大于等于該時(shí)刻液相體積，則將該體單元網(wǎng)格的中心點(diǎn)高度近似看作當(dāng)前時(shí)刻液面 高度，即兩相分界面的位置。如果所有體單元網(wǎng)格尺寸均在毫米量級，則這種方法確定的液 面高度滿足工程設(shè)計(jì)誤差范圍。
[0031] 2.2 S2S輻射模型
[0032] 在本算法中，兩相流區(qū)被節(jié)點(diǎn)化處理，所以在計(jì)算過程中流體網(wǎng)格不會被加載入 求解器，需要使用不受流場影響的輻射模型（步驟104)。表面輻射換熱模型（S2S)在瞬態(tài)計(jì) 算開始之前，需要使用流體網(wǎng)格來計(jì)算輻射角系數(shù)，但在瞬態(tài)計(jì)算中不依賴流場，適合在本 方法中使用。
[0033] 使用S2S模型計(jì)算輻射熱流時(shí)，所需消耗的時(shí)間與參與輻射的面單位數(shù)目的平方 成正比。這意味著，如果以面元為單位計(jì)算熱輻射，需要消耗大量的時(shí)間。所以，可以將若干 個(gè)面元組合成一個(gè)表面束，然后以表面束為單位計(jì)算輻射熱流。
[0034]對于表面束t，其本身的輻射為：
[0035] 為·馬=八·σ_7;4
[0036] 其中，1*為表面束t的輻射率，σ為波爾茲曼常數(shù)，Eb為黑體輻射力，Tt為表面束t的 絕對溫度?？紤]到投射到表面的輻射G會有（1-A t) · Gt的部分被反射，表面束t的有效輻射Jt 可以表示為：
[0037] Jt = At · Eb+(l^t) · Gt
[0038] 講一步的，有有效輻射方稈為：
[0039]
[0040]其中Ft,i為表面束t到表面束i得輻射角系數(shù)。因?yàn)椋?br>[0041]
[0042] 所以有效輻射方程的系數(shù)矩陣嚴(yán)格對角占優(yōu)，故可逆。通過求解該方程，可以獲得 所有表面束的有效輻射熱流密度J。最終，表面束t吸收的輻射熱流密度的表達(dá)式為：
[0043]
[0044]其中Qrt為表面束t包含的所有面元的輻射熱流密度。將輻射熱流密度與通過面元 進(jìn)入固體區(qū)域的對流換熱熱流密度相加，就可以得到通過面元進(jìn)入固體的總熱流密度。
[0045] 由于以表面束為單位計(jì)算輻射角系數(shù)會造成一定誤差，所以在實(shí)際計(jì)算中，需要 在有效輻射方程、輻射熱流密度表達(dá)式的基礎(chǔ)上加以修正。
[0046] 具體包括，