高 亮，王亞軍，邢蘭昌，魏 偉，韓維峰

（1.中國石油大學(xué)（華東）控制科學(xué)與工程學(xué)院，山東青島 266580；2.中國石油勘探開發(fā)研究院新能源研究所，河北廊坊 065007）

0 引言

天然氣水合物在自然界中儲量豐富，具有能量密度高、清潔無污染等優(yōu)點，是一種潛在的戰(zhàn)略性新興能源［1-2］。由于天然氣水合物主要賦存于海底沉積物中，獲取含水合物沉積物實物樣品的成本高且樣品保真度差，因此室內(nèi)水合物模擬實驗成為研究含水合物沉積物物性參數(shù)和儲集參數(shù)的重要手段。實驗中所測量的眾多參數(shù)中，沉積物的水合物飽和度是一個關(guān)鍵的基礎(chǔ)性參數(shù)。目前實驗室內(nèi)水合物飽和度的測量手段主要有電阻率法［3-4］、聲速法［4-5］、CT 掃描技術(shù)［6］、時域反射法（Time Domain Reflectometry，TDR）［7］等。

TDR技術(shù)在土壤物性參數(shù)檢測方面應(yīng)用較多，如對土壤的體積含水量、干密度等進(jìn)行快速測量［8-15］。近年來TDR技術(shù)被用于含水合物沉積物的含水量測量，進(jìn)而實現(xiàn)對水合物飽和度的估算［16］。Zegelin等［17］設(shè)計了三針式TDR探頭，解決了阻抗失配問題，增強了信號的抗干擾能力。Wright 等［7］將三針式TDR探頭用于含水合物沉積物的含水量測量，建立了含水合物沉積物的含水量與表觀介電常數(shù)經(jīng)驗關(guān)系式。Kliner等［18］利用四針式TDR 探頭測量得到了合成制冷劑（R-11）水合物樣本的水合物含量，并建立了介電常數(shù)與水合物含量的關(guān)系式。胡高偉等［19］利用同軸型探頭并在探針上熱敷絕緣套管，實現(xiàn)了對高鹽分模擬海洋沉積物含水量的測量；孫中明［20］利用三針式TDR探頭測量了沉積物中四氫呋喃水合物和甲烷水合物的飽和度。

由電磁波傳播特性可知，垂直于TDR探針長度方向平面上的電場能量在空間上分布不均勻［21］。為了分析不同位置處介質(zhì)對TDR測試響應(yīng)的影響，研究者引入了空間權(quán)重函數(shù)的概念，并用來表示TDR探測的空間靈敏度［22-25］。含天然氣水合物沉積物的介電常數(shù)在空間中也呈現(xiàn)出顯著的不均勻分布的特性［26］。由于TDR探頭具有特殊的幾何形狀以及被測介質(zhì)（如含水合物沉積物）在空間分布呈現(xiàn)非均勻性，因此用解析方法分析被測介質(zhì)對TDR 響應(yīng)的影響規(guī)律存在較大的局限性。本文對垂直于TDR 探針方向平面上測試區(qū)域內(nèi)含水合物沉積物的測量響應(yīng)進(jìn)行了有限元數(shù)值建模和分析，基于數(shù)值模型對TDR探頭的結(jié)構(gòu)參數(shù)進(jìn)行了優(yōu)化，利用優(yōu)化后的探頭對水合物飽和度進(jìn)行了仿真測量，并與表觀介電常數(shù)與水合物飽和度關(guān)系模型進(jìn)行了比較，從而為用于含水合物沉積物的TDR探頭的結(jié)構(gòu)參數(shù)優(yōu)化和水合物飽和度反演模型的建立提供模型基礎(chǔ)和方法參考。

1 TDR原理與數(shù)值建模分析

1.1 工作原理

TDR技術(shù)利用電磁波在介質(zhì)中的傳播速度來確定被測介質(zhì)的表觀介電常數(shù)Ka。在被測介質(zhì)中，電磁波傳播速度vp與Ka的關(guān)系為

式中，c為真空中電磁波傳播速度，3.0 ×108m／s。典型的三針式TDR探頭結(jié)構(gòu)如圖1 所示，探針長度為L，假設(shè)電磁波抵達(dá)探針首端時刻為t1，抵達(dá)探針末端后返回至探針首端時刻為t2，則t1與t2的時間差Δt為電磁波沿TDR探頭傳播的總時間，由此可得電磁波在被測介質(zhì)中的傳播速度為

圖1 三針式TDR探頭示意圖

將式（1）、（2）聯(lián)立可得［27］：

1.2 數(shù)值建模

（1）幾何結(jié)構(gòu)與材料。忽略電磁波沿TDR 探頭傳播方向的衰減以及探頭的末端效應(yīng)，針對垂直探針長度方向的截面進(jìn)行研究。將該截面上的研究區(qū)域設(shè)置為半徑為R的圓（見圖2），內(nèi)探針設(shè)置在圓心處，外探針中心與內(nèi)探針中心的距離均為s，內(nèi)／外探針的直徑均為d。外探針的位置根據(jù)探針的個數(shù)進(jìn)行設(shè)定，圖2 所示為探針個數(shù)分別為3、4 和5 的情形。

圖2 三／四／五針式TDR探頭的探針平面分布

針對被測介質(zhì)為含水合物沉積物的情形，在上述研究區(qū)域內(nèi)模擬由20～30 目海砂構(gòu)成的多孔介質(zhì)。多孔介質(zhì)的幾何結(jié)構(gòu)和對計算區(qū)域進(jìn)行剖分得到的網(wǎng)格如圖3 所示，其中圓形區(qū)域代表海砂顆粒，四邊形區(qū)域代表水合物顆粒，探針為三針式TDR探頭中的內(nèi)探針，剩余區(qū)域設(shè)置為水。模型的網(wǎng)格形狀為自由三角形，由物理場控制網(wǎng)格類型。模型中孔隙水、海砂顆粒和水合物的相對介電常數(shù)分別為εw＝80、εs＝4.0 和εh＝7.0［28］。

圖3 含水合物沉積物模型的幾何結(jié)構(gòu)及其網(wǎng)格

（2）控制方程與邊界條件。在垂直于探針長度方向的截面內(nèi)，TDR探針周圍的電磁場滿足電場高斯定律和本構(gòu)方程且電荷密度為0：

式中：D 為電通量強度（C／m2）；ρ 為電荷密度（C／m3），此處取值為0；E為電場強度（V／m）；ε0和εr分別為真空介電常數(shù)和相對介電常數(shù)。

聯(lián)立此二方程得：

由于電場強度E為電勢Φ的梯度，即

式中，Φ 是標(biāo)量函數(shù)，表示垂直于探針長度方向的截面內(nèi)的電勢分布，滿足拉普拉斯方程［29］。

將式（7）代入式（6），得

通過求解以上方程，可得到截面內(nèi)各空間位置的電勢和電場強度。

模型采用有限元計算軟件COMSOL 中AC／DC 模塊中“電場和電流”下的“靜電”物理場。如圖4 所示，3 個探針分別代表3 個終端，內(nèi)探針終端所在邊界設(shè)置為恒定電勢1 V，兩根外探針終端所在邊界設(shè)置為恒定電壓勢-1 V，研究區(qū)域的外邊界為探針長度方向垂直截面的外邊界（見圖4），設(shè)置為零電荷。

圖4 數(shù)值模型中三針式TDR探頭與研究區(qū)域

（3）模型驗證。首先將模型中被測介質(zhì)部分全部設(shè)置為空氣，垂直于探針截面上的電勢分布如圖5 所示。由圖可見：內(nèi)探針周圍電勢等高線由內(nèi)至外遞減至-1 V，外探針周圍的電勢等高線由外至內(nèi)遞減至-1 V。該等電勢分布結(jié)果與Zegelin［17］和詹良通［25］等通過求解拉普拉斯方程繪制的探針之間的電勢線計算結(jié)果一致，從而證明本文所建立數(shù)值模型正確。

圖5 垂直于探針截面上電勢和電場強度分布

1.3 數(shù)據(jù)分析方法

首先通過數(shù)值模型獲得研究區(qū)域的電勢分布，然后利用電勢分布計算空間權(quán)重函數(shù)，最后計算研究區(qū)域被測介質(zhì)的表觀介電常數(shù)，即為TDR測得介電常數(shù)的仿真計算值。垂直于TDR 探針方向截面上的電勢分布Φ（x，y）由該截面上介電常數(shù)的分布所決定，其中（x，y）表示界面上任意點的位置。TDR探頭周圍每個空間坐標(biāo)點的空間權(quán)重函數(shù)為［30］：

式中：Ω代表研究區(qū)域；Φ0（x，y）是研究區(qū)域中具有均勻分布介電常數(shù)被測介質(zhì)時的電勢分布；Φ（x，y）是研究區(qū)域中為實際被測介質(zhì)時的電勢分布?？臻g權(quán)重函數(shù)描述了TDR探頭在研究區(qū)域Ω內(nèi)任何點的測量靈敏度?？臻g權(quán)重函數(shù)滿足：

由式（9）計算得到的表觀介電常數(shù)為

式中，K（x，y）是研究區(qū)域Ω中點（x，y）處被測介質(zhì)的相對介電常數(shù)。

利用研究區(qū)域的電勢分布可以計算出探針周圍電場的有效極化面積。有效極化面積是指對TDR 測量響應(yīng)有主要貢獻(xiàn)區(qū)域的面積，在該區(qū)域之外，被測介質(zhì)性質(zhì)的改變對TDR測量響應(yīng)沒有顯著的影響［24］?；诳臻g權(quán)重函數(shù)的概念，定義垂直探針方向截面上的部分區(qū)域?qū)DR 測量響應(yīng)的貢獻(xiàn)在總響應(yīng)中所占的比例f，即

在數(shù)值模型中，式（12）中分子為部分有限元網(wǎng)格區(qū)域內(nèi)網(wǎng)格的空間權(quán)重函數(shù)wi與網(wǎng)格面積Ai乘積的和。從內(nèi)探針周圍開始，按照空間權(quán)重函數(shù)從大到小的順序進(jìn)行排序，whi為研究區(qū)域內(nèi)空間權(quán)重函數(shù)的最大值。當(dāng)f＝90%時，認(rèn)為此時區(qū)域面積為TDR的有效極化面積，當(dāng)f＞90%時，區(qū)域面積的增加對TDR測量響應(yīng)的影響不顯著。

2 TDR探頭參數(shù)優(yōu)化

2.1 探針的個數(shù)

為探究探針的個數(shù)對TDR測量響應(yīng)的影響，仿真分析了水合物位于測量區(qū)域中不同位置時各探針周圍的電場強度分布以及表觀介電常數(shù)的變化。

針對具有不同探針個數(shù)的TDR探頭（探針個數(shù)分別為3、4、和5），建立了垂直于探針長度方向截面的多孔介質(zhì)二維數(shù)值模型。模型的幾何結(jié)構(gòu)如圖2 和3 所示，相關(guān)幾何參數(shù)取值為R＝20 mm、s＝15 mm、d＝3 mm。此時二維模型的總面積為1 256 mm2，沙粒和孔隙的面積分別為779 mm2和477 mm2，孔隙度φ為0.38。圖6 為在探針周圍3 個典型區(qū)域分別加入占據(jù)相同面積的水合物，從而模擬含水合物多孔介質(zhì)，水合物區(qū)域的面積為30.5 mm2，占模型總面積的2.43%。以三針式TDR 探頭為例，圖6 中橙色圓點代表海砂顆粒，區(qū)域A在內(nèi)外探針之間，區(qū)域B在外探針外側(cè)，區(qū)域C在內(nèi)探針上方，四針式、五針式探針模型的水合物區(qū)域位置與三針式探針模型一致，見圖7 中的方框。

圖6 水合物在三探針TDR探頭模型中的空間位置

對于具有不同探針個數(shù)的探頭，以水合物位于區(qū)域A為例繪制出探針周圍的電場強度分布，如圖7 所示。由于TDR探頭的內(nèi)外探針間距相同，由電場強度分布圖可知，隨著探針個數(shù)由3 增加至5，高電場覆蓋區(qū)域形狀由不規(guī)則橢圓形變?yōu)橐酝馓结槥轫旤c的三角形和四邊形，有效極化面積逐步增大。

圖7 不同探針個數(shù)時水合物位于區(qū)域A時電場強度

水合物位于不同區(qū)域時測得表觀介電常數(shù)的相對變化量（相對于無水合物的情況）與探針個數(shù)的關(guān)系如圖8 所示。對于不同探針個數(shù)的探頭，當(dāng)水合物位于區(qū)域A時，Ka相對變化量均保持最大，可見位于內(nèi)外探針之間的水合物對測量的Ka影響最大，而位于外探針之外的水合物對Ka影響最小。隨著探針個數(shù)的增加，水合物位于區(qū)域A和B時測得Ka相對變化量逐漸減小，位于區(qū)域C時測得Ka相對變化量逐漸增大?？傮w上，三探針TDR 探頭測量得到的Ka變化量較大，即測量靈敏度相對更高。在實際應(yīng)用中，當(dāng)較多的探針插入被測沉積物體系，其將對被測體系結(jié)構(gòu)產(chǎn)生較大的擾動，不利于對原狀樣品進(jìn)行參數(shù)測量。綜合考慮，本文采用三探針式TDR探頭。

圖8 水合物位于不同區(qū)域時測得表觀介電常數(shù)的相對變化量與探針個數(shù)的關(guān)系

2.2 探針的間隔與直徑

電場能量密度正比于電場強度的平方，TDR 探針周圍的電場能量密度越大，則該位置處的介質(zhì)對測試結(jié)果的影響越大，利用探針之間電場能量密度的分布可以分析TDR的測量響應(yīng)。探針間距和探針直徑對電場能量密度分布的影響并不相互獨立，當(dāng)內(nèi)外探針之間距離過大或者探針本身直徑過小時，電場能量將被限制在探針附近的區(qū)域，即出現(xiàn)“集膚效應(yīng)”，此時電場有效極化面積將大幅縮小，導(dǎo)致有效探測區(qū)域面積的減小。為了緩解“集膚效應(yīng)”所導(dǎo)致的電場能量在探針附近過度集中的問題，Knight 建議探針直徑應(yīng)大于探針間距的十分之一［22］。

假設(shè)圓筒型的水合物反應(yīng)釜橫截面直徑為40 mm，參照圖2 將內(nèi)外探針之間的距離和探針直徑分別表示為s和d。已有研究表明，如圖7（a）中，三針式TDR探頭在垂直于探針長度方向截面內(nèi)形成橢圓形電場，橢圓長軸的一半比s+d／2 大3 mm 左右［25］，并且d應(yīng)大于s／10［22］，本文首先設(shè)定內(nèi)外探針之間距離為15 mm。在確定s的前提下，為了探究探針直徑對TDR測量響應(yīng)的影響，仿真分析了d從s／10 逐漸增大時，探針周圍最大電場能量密度和有效極化面積的變化規(guī)律。以下建立了一系列TDR探頭模型，探針直徑從1.5 mm開始，依次增加0.5 mm。

圖9 所示為3 種d的條件下探針周圍的電場能量密度分布，可見電場能量密度最大值出現(xiàn)在內(nèi)探針邊緣和外探針靠近內(nèi)探針一側(cè)的邊緣。最大能量密度越大則顯示越強的“集膚效應(yīng)”，即電場能量更多地被限制于探針附近的區(qū)域。圖10 展示了最大電場能量密度和有效極化面積隨探針直徑變化的規(guī)律。由圖10可知，隨著d的逐漸增大，最大電場能量密度逐漸減小、有效極化面積逐漸增大。最大電場能量密度值越小，則“集膚效應(yīng)”越弱，有效極化面積越大，探針的探測范圍越大；但是探針直徑越大，探針橫截面積占據(jù)兩針之間面積的比例越大，則對被測樣品產(chǎn)生的干擾越大［31］。分析圖10 曲線的轉(zhuǎn)折點位于d＝3 mm附近，本文確定本研究的三針式探頭的內(nèi)外探針距離為15 mm、探針直徑為3 mm。

圖9 不同直徑時探針電場能量密度分布

圖10 不同探針直徑時探針周圍最大電場能量密度和有效極化面積

3 水合物飽和度仿真測量

通過改變多孔介質(zhì)數(shù)值模型中水合物區(qū)域的大?。ㄒ妶D3），分別模擬0%～100%水合物飽和度下的含水合物沉積物，并仿真分析水合物飽和度與表觀介電常數(shù)的關(guān)系。其中，表觀介電常數(shù)用式（11）計算得到。

以數(shù)值模型的計算結(jié)果為基礎(chǔ)，對水合物飽和度Sh與Ka進(jìn)行擬合得到：

將Sh與體積含水量θV之間的關(guān)系式［16］（即：Sh＝1 -θV／φ）代入式（13）可得

圖11 比較了本文擬合關(guān)系式（14）與Topp 經(jīng)驗式［10］、Wright 經(jīng)驗式［7］、CRIM（Complex Refractive Index Method）模型［32］和MD（Maxwell-DeLoor）模型［33］。

Topp經(jīng)驗式：

Wright經(jīng)驗式：

CRIM模型：

MD模型：

分析圖11 可知：①基于數(shù)值仿真數(shù)據(jù)擬合得到的關(guān)系式（14）、MD模型與仿真數(shù)據(jù)一致性較高，在0～100%水合物飽和度范圍內(nèi)兩者的均方根誤差分別為2.92 和3.74%；②CRIM模型適用于水合物飽和度較高的情況，如Sh＞50%，在該范圍內(nèi)CRIM 模型的均方根誤差為2.49%；③Topp 經(jīng)驗式和Wright 經(jīng)驗式對仿真數(shù)據(jù)的適應(yīng)性較差，兩者均方根誤差分別為12.52 和10.57%。

圖11 本文仿真擬合與其他經(jīng)驗式或模型對比圖

4 結(jié)語

針對TDR探頭在含水合物沉積物中的應(yīng)用，基于有限元數(shù)值模型計算了垂直于探針方向平面上測試區(qū)域內(nèi)含水合物沉積物的TDR 響應(yīng)，在對探針個數(shù)、探針間距和直徑等結(jié)構(gòu)參數(shù)進(jìn)行優(yōu)化設(shè)計的基礎(chǔ)上，對多孔介質(zhì)中水合物飽和度進(jìn)行了仿真測量。得到以下三點認(rèn)識：

（1）隨著探針個數(shù)由3 增加至5，高電場覆蓋區(qū)域形狀由不規(guī)則橢圓形變?yōu)橐酝馓结槥轫旤c的三角形和四邊形，有效極化面積逐步增大；三針式TDR 探頭對于內(nèi)外探針之間的水合物探測靈敏度最高。

（2）在探針間距一定的條件下，隨著探針直徑的增大，最大電場能量密度減小，表明“集膚效應(yīng)”越弱，有效極化面積增大，則探測范圍越大；但探針直徑越大，其對被測樣品產(chǎn)生的干擾越大。

（3）基于有限元數(shù)值模型得到的含水合物沉積物TDR響應(yīng)與Maxwell-DeLoor 模型具有較好的一致性，驗證了三針式TDR 探頭測量響應(yīng)數(shù)值仿真計算的正確性。依據(jù)本文所優(yōu)化的探針結(jié)構(gòu)參數(shù)加工制作TDR探頭，進(jìn)一步開展含天然氣水合物模擬沉積物的TDR實驗測試，可以獲得不同水合物飽和度條件下的表觀介電常數(shù)實測數(shù)據(jù)，進(jìn)而綜合利用有限元數(shù)值模擬和實驗測試數(shù)據(jù)來建立可靠性高、準(zhǔn)確度高、適用范圍廣的水合物飽和度計算模型。