基于LTNE的飽和多孔介質THMC耦合模型數值研究

戴 倩， 李培超

(上海工程技術大學 機械與汽車工程學院， 上海 201620)

無論是生產還是生活中，多孔介質材料都十分常見，例如海綿、隔音泡沫及鋰離子電池隔膜等。多孔介質是由固體骨架和孔隙空間組合構成的物質，孔隙一般由單相或多相流體共同填充。其研究和發展擁有較為悠久的歷史[1]。為簡單起見，前人采用了一些通俗而經典的物理模型進行研究，如圓環[2]、圓管[3-4]、平板[5]和無限大平面[6]等，以軟件輔助模擬，用公式推導分析，大大促進了多孔介質理論的進步。

對于多孔介質的研究，大多集中于多物理場(諸如應力應變場、流場、溫度場、電化學場等)耦合問題中。常見的有流固耦合[7-8]、熱流耦合[9]和熱流固耦合[10-13]等數學模型。但更完善的為熱流固化耦合模型，常應用于生物醫療、海水入侵、水利工程和農田噴灑等領域。

諸多科研工作者對多孔介質熱流固化耦合數學模型進行了深入探討和研究。孔祥言[14]404建立了天然氣水合物的部分耦合控制方程組，并考慮了滲透率和毛管力隨飽和度的變化。Sun等[15]描述了熱流固化耦合力學行為，將模擬結果與測試數據進行了對比和分析，驗證新開發代碼的實用性。盛金昌等[16]240著重考慮鉆井過程中熱流固化耦合過程和井壁孔隙壓力、溫度、應力和離子溶質摩爾分數等變化，提供了更加精確的井壁應力分布。

研究多孔介質多物理場耦合問題的傳統假設大都是基于局部熱平衡理論即LTE[17-18]，但當內部換熱不均勻，固體骨架和孔隙流體的溫度差值不可忽略時，局部熱平衡理論不再適用，需運用更加合理的局部熱非平衡理論即LNTE[19-20]來描述此類情形。例如，陳麗等[21]750基于LTNE理論和加權平均溫度的概念，把握化學-熱-孔隙彈性模型的特點，應用Laplace變換法對多場耦合表達式進行無量綱化，并與LTE條件的相關結果進行圖像對比，發現LTNE效應在Biot數為中等范圍內不可忽略。

前文所述的多物理場耦合方式，多為部分耦合或弱耦合。鑒于此，課題組假設多孔介質中的孔隙流體為溶質(NaCl)和溶劑(H2O)2種化學組分，構建更加完備的THMC強耦合數學模型。

1 數學模型

課題組將空間問題簡化為平面問題，如圖1所示。采用內含一半徑a=0.1 m的圓孔無限大平面飽和多孔介質模型進行研究。其中，格子區域為飽和多孔介質無限大平面，表示以r的圓心為原點的徑向坐標。簡單起見，基本假設如下：①多孔介質為均勻且各項同性的介質，固體骨架變形符合小變形假設；②流體為單一的、穩定的和不可壓縮的流體，流體流動遵從Darcy定律；③忽略自然對流、擴散和輻射換熱的影響；④拉應力為正。

圖1 內含圓孔的多孔介質無限大平面示意圖Figure 1 Schematic diagram of porous medium infinite plane containing circular hole

1.1 控制方程

基于孔祥言[14]《高等滲流力學》第3版(第1，7，12章)，李培超等[7]422飽和多孔介質流固耦合滲流的數學模型，陳麗等[21]752化學-熱-彈性模型和本課題組里岳飛龍飽和多孔介質單相熱流固完全耦合模型，建立基于LTNE理論的THMC完全耦合的控制方程組。

為降低模型的復雜程度，加快軟件的求解速率，將控制方程組簡化為一維軸對稱的情形進行求解，即將笛卡爾坐標轉換為極坐標形式(自變量x，y和z轉換成自變量r)。

因此，關于一維軸對稱坐標下的熱流固化完全耦合控制方程為：

(1)

(2)

(3)

(4)

(5)

式中：ε為多孔彈性體的應變，p為孔隙流體壓力，m為溶質摩爾分數變化量，λ為Lame常數，G為Lame常數，K為體積彈性模量，γ為化力耦合參數，β為體熱膨脹系數，α為Biot系數，Ts為固體骨架溫度，Tf為孔隙流體溫度，h為流固界面換熱系數，kTs為固體導熱率，kTf為流體導熱率。每一個物理場控制方程都含有一個主要因變量。

1.2 初始條件和邊界條件

數學模型需要一定的定解條件方可求解，即需建立初始條件和邊界條件。

設定各物理場在a≤r<∞區域范圍內的初始條件為：

(6)

各物理場的邊界條件如下：

(7)

此處假設，ma和Ta分別為圓孔邊界r=a處溶質摩爾分數的變化增量和固體與流體的加載溫度，ma=1.8×10-2，Ta=50 ℃，H(t)為Heaviside函數。

其中，根據文獻[22]第12頁的內容，整理推導出相關的耦合系數和參數方程如下：

(8)

2 數值方法

2.1 COMSOL Multiphysics求解

COMSOL內部已具有較為完善的模塊，但功能最強大的、運用最靈活的仍是其PDE(偏微分方程)模塊[16]242。PDE模塊可求解用戶自定義所研究問題的復雜控制方程組，根據系統提供的假設方程輸入，進行數值模擬。

2.2 參數取值

表1為模型中各物性參數取值，部分參數取自參考文獻[16]～[17]和參考文獻[21]～[22]。

表1 固體和流體的物性參數

3 結果與討論

3.1 驗證

課題組將熱流固化4場耦合退化為熱流固化3場耦合，通過與陳麗等[21]758的數據對比，驗證本研究數學模型的可行性和數值結果的可靠性，如圖2所示。圖中數據一致性一定程度上說明本方案的可實施性。

圖2 本研究退化模型與參考文獻模型對比Figure 2 Comparison of proposed degradation model and reference model

3.2 結果分析

課題組借助COMSOL軟件平臺得到多物理場耦合數值模擬結果。為更好呈現結果，課題組將因變量和自變量進行無量綱化處理。

(9)

式中：mD為量綱為一的摩爾分數變化量，pD為量綱為一的孔隙壓力，εD為量綱為一的應變，θs為量綱為一的固體溫度。

圖3為化學場隨時空變化分布圖。溶質摩爾分數變化量分布不明顯，課題組截取部分數據放大研究(下文有關溶質摩爾質量分數變化量分布均采取局部放大圖進行研究)。由圖(b)可明顯看出，溫度升高會引起溶質微妙變化，故剛開始的邊界處，溶質摩爾分數變化量有所波動，但波及范圍較小，后隨時間遞增和溶質擴散，波及范圍逐漸擴張。

圖3 化學場隨時空變化分布圖Figure 3 Diagram of chemical fields varying over time and space

3.3 參數分析

由于本研究中控制方程中參數較多，故選取有代表性的參數進行研究，分析其對結果影響深淺。課題組選取2 000 s時，研究滲透率k，固體導熱率kTs和Biot系數α對多物理場耦合作用。

3.3.1 滲透率的影響

滲透率表征流體在多孔介質中滲透能力的大小。圖4展示了不同滲透率對結果的影響，從圖中可以看出，孔隙壓力隨著滲透率的增大先減小后增大。這是由于滲透能力越強，流體流動性越好，孔隙壓力能夠得到更快地釋放，故會出現滲透率越大孔隙壓力越小的現象，但到達一定位置后，孔隙壓力會回歸之前狀態。

圖4 t=2 000 s時滲透率對孔隙壓力影響Figure 4 Influence of permeability on pore pressure at t=2 000 s

3.3.2 固體導熱率的影響

導熱率是材料最重要的熱物性參數之一，描述其導熱性能。其對化學場和溫度場的影響，如圖5所示。由圖可知，隨固體導熱率增加，固體溫度增加，溫度變化會使得溶質摩爾分數發生細微的變化。

圖5 t=2 000 s時固體導熱率對mD和θs的影響Figure 5 Influence of heat transfer rate to mD and θs at t=2 000 s

3.3.3 Biot系數的影響

Biot系數又稱有效應力系數，描述因孔隙流體支撐載荷而引發骨架剛度變化程度。從圖6可看出Biot系數對前幾個物理場作用較明顯，這是由于骨架剛度變化會引起應力場發生變化，從而引發其他物理場的變化，但應變對溫度場貢獻較小，Biot系數對溫度場影響微乎其微，故此處未展示溫度場的分布。

圖6 t=2 000 s時Biot系數的影響Figure 6 Influence of Biot coefficient at t=2 000 s

4 結論

課題組基于LTNE理論在多孔介質內發生THMC耦合行為的影響，得出結論如下：

1) 在設定的初邊值條件下，隨著時間的增加，各物理場相互作用和影響，變化范圍逐漸擴大；尤其是化學場，其摩爾分數的影響由深到淺，但最終都將趨于穩定。

2) 各參數對THMC耦合模型影響深遠，需選取適當參數進行研究(如對于參數選值，考慮是否會產生較高的熱應力，防止結構破壞)；與其他參數比較，導熱率對溫度的影響較大。

3) 圓孔邊界0.1～0.2 m處THMC耦合作用十分明顯，說明流體在孔隙中釋放的反應速率較快，并周向傳播；且在此處，溫度場對各物理場的貢獻較顯著，但反之影響較小。