戴 倩， 李培超

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

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

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

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

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

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

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

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

圖1 內(nèi)含圓孔的多孔介質(zhì)無(wú)限大平面示意圖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ù)雜程度，加快軟件的求解速率，將控制方程組簡(jiǎn)化為一維軸對(duì)稱的情形進(jìn)行求解，即將笛卡爾坐標(biāo)轉(zhuǎn)換為極坐標(biāo)形式(自變量x，y和z轉(zhuǎn)換成自變量r)。

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

(1)

(2)

(3)

(4)

(5)

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

1.2 初始條件和邊界條件

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

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

(6)

各物理場(chǎng)的邊界條件如下：

(7)

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

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

(8)

2 數(shù)值方法

2.1 COMSOL Multiphysics求解

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

2.2 參數(shù)取值

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

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

3 結(jié)果與討論

3.1 驗(yàn)證

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

圖2 本研究退化模型與參考文獻(xiàn)模型對(duì)比Figure 2 Comparison of proposed degradation model and reference model

3.2 結(jié)果分析

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

(9)

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

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

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

3.3 參數(shù)分析

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

3.3.1 滲透率的影響

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

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

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

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

圖5 t=2 000 s時(shí)固體導(dǎo)熱率對(duì)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ù)對(duì)前幾個(gè)物理場(chǎng)作用較明顯，這是由于骨架剛度變化會(huì)引起應(yīng)力場(chǎng)發(fā)生變化，從而引發(fā)其他物理場(chǎng)的變化，但應(yīng)變對(duì)溫度場(chǎng)貢獻(xiàn)較小，Biot系數(shù)對(duì)溫度場(chǎng)影響微乎其微，故此處未展示溫度場(chǎng)的分布。

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

4 結(jié)論

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

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

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

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