多孔介質中單相氣體局部流動的均質化建模

李樹光, 曲 凱

(大連海事大學 理學院, 遼寧 大連 116026)

0 引 言

多孔介質中的氣體流動廣泛存在于自然現象和工業技術中,因而對其開展數學建模研究具有重要的科學意義和應用價值[1],如地質層中天然氣的輸運[2],飛機油箱內多孔閥的油氣過濾,飛行器外殼熱防護材料的氣動燒蝕等．然而,由于測量技術的局限性和多孔介質空間結構的復雜性,使得這種流動的研究非常困難,大多數工作僅能獲得多孔介質中氣體流動的宏觀結果．許多研究是根據Darcy定律及其修正,從宏觀現象分析孔隙中氣體的流動．這種情況下,單孔中所發生的局部效應僅作為整體特征的簡單估計(如Darcy定律中平均流速和實驗中滲透系數的確定),未充分考慮孔隙中流動的真實過程．孔隙中的流動對滲流的整體影響非常重要[3]:孔隙的幾何特征和局部流動決定了滲流的滲透性．因此,有必要準確描述孔隙中的局部流動．

許多學者應用多種方法研究了多孔介質孔隙中的局部流動,如均質化和體積平均模型[4],孔隙網絡(“管束”)模型[5],分數階模型[6]等．每種方法都有其優缺點,對于滲流問題的研究,描述孔隙中發生的實際物理過程以及微宏觀過程之間的理論聯系至關重要．因此,本文采用了文獻[7]中提出的針對周期性復合材料建模的漸近均質法,該方法能夠有效模擬多孔介質孔隙中的物理過程．目前,漸近均質法已經被廣泛應用于解決多孔材料的熱力學問題[8]、復合材料有效性能預測[9-10],以及材料滲透性能估計[11]等問題．

應用漸近均質法獲得的描述局部流動的數學模型具有一些特殊性[12],主要是:通常具有積分-微分類型的均勻化限制和周期性邊界條件所構成的約束條件．這使得局部問題的理論意義和求解方法有別于經典的數學模型．一些研究數值求解了周期單元上的局部問題,包括罰函數有限元法[13]、光滑粒子流體動力學(SPH)方法[14]等．但是少有討論局部問題的特殊性和復雜約束條件的求解方法．

本文的目的是研究周期單元上氣體的局部流動問題,提出了一種基于對稱性和反對稱性擴展的簡化方法,開發求解局部問題的最小二乘有限元法,通過微管中Poiseuille流動的解析解驗證了相關模型和算法．

1 多孔介質中氣體流動的連續介質力學模型

1.1 多孔介質的幾何模型

考慮如圖1所示的正交各向異性多孔介質．建立如下假設: ① 多孔介質的微觀結構具有周期性特征; ② 多孔介質無死孔; ③ 周期單元是孔隙中充滿氣體的纖維和纖維材料所組成的區域; ④ 周期單元關于局部直角坐標系中的坐標平面具有幾何和物理的對稱性,這使得能夠將計算區域簡化為1/8周期單元．

(a) 多孔介質的微觀結構 (b) 多孔介質的周期單元 (c) 多孔介質的1/8周期單元 (a) The microstructure of the porous medium (b) Periodic cell of the porous medium (c) 1/8 periodic cell of the porous medium

圖1中,V表示整個孔隙所占的區域,其邊界(固-氣邊界)用Σ表示,Vξ表示周期單元,Σξ表示孔隙中氣體與固體的邊界．

1.2 氣體流動的力學方程組

多孔介質中氣體的流動由可壓縮Navier-Stokes方程組描述:

(1)

其中,ρ是密度,v是速度,T是Cauchy應力張量,?表示張量積．

對氣體介質設置如下假設: ① 孔隙中的氣體介質為理想氣體; ② 氣體是各向同性的Newton黏性介質,黏性系數非常小且不為零; ③ 氣體的質量力密度為零; ④ 氣體的運動是等溫等壓的; ⑤ 氣體的體積黏度系數[12]為零,即λ+2μ/3=0,其中λ和μ是黏度系數．

根據各向同性線黏性介質的本構關系,Cauchy應力張量可寫為

T=-pE+σ,

(2)

其中,p是壓力,σ=λ(?·v)E+2με是黏性應力張量,ε=(??v+??vT)/2是應變率張量．

考慮理想氣體的狀態方程:

p=ρRθ0,

(3)

其中,R是氣體常數;θ0是恒定溫度．

在固氣界面,設置無滑移邊界條件:

v|Σ=0．

(4)

初始時刻t=t0的給定壓力為

p|t=t0=p0．

(5)

2 多孔介質中氣體流動的局部建模

2.1 局部建模中漸近均質法的尺度放大

根據多孔介質的周期性假設,速度v、壓力p和密度ρ由小參數κ的漸近展開式表示:

(6)

根據氣體黏性的假設,可以得到μ=μ0κ2．

將式(6)和黏性假設代入系統(1)—(5)并以小參數κ的冪次進行分組后,得到周期單元上的局部問題如下:

(7)

進一步考慮三維多孔結構,氣體沿著局部坐標軸運動時,由于局部問題(7)的線性性質,其解可表示為“輸入數據”?xp(0)的線性函數:

(8)

(9)

(10)

2.2 局部問題的物理意義

通過分析系統(9)可得結論:在局部層面上,當系統(9)固定一個α的值,也就確定一個偽不可壓縮線性流體的穩態流動問題,并且系統(9)的解僅由孔的幾何形狀決定．因此,系統(9)適用于在上述假設框架下計算多孔介質中氣體的流動．

2.3 局部問題在1/8周期單元上的簡化

進一步,考慮將局部問題簡化為1/8周期單元上的問題,從而進行數值求解．因此,根據解的對稱和反對稱原理得到以下關于解的延續性定理[11-12]．

(11)

2,3時,可采用類似的置換方式進行．圖2中的符號“+”表示從第一象限擴展到其他象限時函數不更改符號(即對稱擴展),符號“-”表示符號更改(即反對稱擴展)．

利用系統(11)的解進行對稱和反對稱擴展,在1/8周期單元的邊界平面上設置如下邊界條件:

(12)

2.4 滲透性的確定

(13)

3 局部問題的最小二乘有限元數值解法

(14)

(15)

(16)

接下來,考慮有限元公式(14)—(16)的近似函數．由于最小二乘公式(14)中近似函數的最小要求是線性的,因此可采用C0基函數[17]．盡管在所用的近似空間之間沒有兼容性限制(即inf-sup條件),在本文中仍然使用滿足inf-sup條件的近似函數．

4 數 值 結 果

本節中,首先驗證了所提出的模型和算法的可靠性和準確性,然后模擬了在圖1所示結構中的流動．

4.1 微管中Poiseuille流動的數學模型和數值模擬

為驗證模型和算法,考慮單通道結構中α=1時模型(9)的解析解．微管的1/8周期單元如圖3所示．

圖3 圓柱形單通道結構的1/8周期單元Fig. 3 The 1/8 periodic cell of the cylindrical single channel structure

(17)

(18)

圓柱形單通道結構的1/8周期單元中,系統(18)在柱坐標系下可表述為

(19)

其中,R是圓柱形單通道的半徑．根據式(19)可得Cauchy問題:

(20)

在試驗中,微管的半徑為0.4．因此,孔隙度的精確結果為0.502 654 82．計算微管中Poiseuille流動采用的是2 674個單元、4 698個節點的四面體網格．圖4中給出了計算結果,孔隙度的數值結果為0.502 654 65．由此,數值結果驗證了模型和數值算法的可靠性和準確性．

(a) Poiseuille流動的數值解 (b) 與精確解的比較 (a) The numerical solution of the Poiseuille flow (b) Comparison of numerical and exact solutions

4.2 正交各向異性多孔介質中局部流動的數值結果

考慮圖1所示的幾何模型,孔隙中球體的無量綱半徑為0.35,圓柱(連接通道)的無量綱半徑為0.1．模型采用的參數為:壓力量綱p0=106Pa,速度量綱v0=1 m/s,多孔介質的量綱長度L=1 m,周期單元的量綱長度l=10-4m,氣體(氮氣)的黏度μ=0.018 448 Pa·s．

表1 多孔介質孔隙率和滲透系數的計算結果

(a) 速度分量 速度分量

5 結 論

本文應用漸近均質法建立了描述三維周期性多孔結構孔隙尺度下單相氣體局部流動的數學模型,獲得了周期單元上的局部問題,明確了局部問題的數學特殊性和物理意義,提出了一種基于對稱性和反對稱性擴展的簡化方法．通過對局部問題的分析,能夠準確地獲得多孔介質的滲透率．

本文在分析局部問題的基礎上,提出了求解局部問題的最小二乘有限元方法,克服了由平均算子和周期性邊界條件引起的數值困難．最后,通過理論分析獲得了微管中Poiseuille流動的解析解,并驗證了所提出的數學模型和數值算法的可靠性和準確性．研究結果表明,應用漸近均質法能夠解決多孔介質中的單相氣體流動問題,獲得多孔介質的氣體滲透性．