基于LTNE的飽和多孔介質(zhì)THMC耦合模型數(shù)值研究

戴 倩， 李培超

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

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

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

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

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

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

1 數(shù)學(xué)模型

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

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

1.1 控制方程

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

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

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

(1)

(2)

(3)

(4)

(5)

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

1.2 初始條件和邊界條件

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

設(shè)定各物理場在a≤r<∞區(qū)域范圍內(nèi)的初始條件為：

(6)

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

(7)

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

其中，根據(jù)文獻[22]第12頁的內(nèi)容，整理推導(dǎo)出相關(guān)的耦合系數(shù)和參數(shù)方程如下：

(8)

2 數(shù)值方法

2.1 COMSOL Multiphysics求解

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

2.2 參數(shù)取值

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

表1 固體和流體的物性參數(shù)

3 結(jié)果與討論

3.1 驗證

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

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

3.2 結(jié)果分析

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

(9)

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

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

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

3.3 參數(shù)分析

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

3.3.1 滲透率的影響

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

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

3.3.2 固體導(dǎo)熱率的影響

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

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

3.3.3 Biot系數(shù)的影響

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

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

4 結(jié)論

課題組基于LTNE理論在多孔介質(zhì)內(nèi)發(fā)生THMC耦合行為的影響，得出結(jié)論如下：

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

2) 各參數(shù)對THMC耦合模型影響深遠，需選取適當(dāng)參數(shù)進行研究(如對于參數(shù)選值，考慮是否會產(chǎn)生較高的熱應(yīng)力，防止結(jié)構(gòu)破壞)；與其他參數(shù)比較，導(dǎo)熱率對溫度的影響較大。

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