基于有限元模型的TDR探頭優化及水合物飽和度測量的數值實驗

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

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

0 引言

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

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

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

1 TDR原理與數值建模分析

1.1 工作原理

TDR技術利用電磁波在介質中的傳播速度來確定被測介質的表觀介電常數Ka。在被測介質中，電磁波傳播速度vp與Ka的關系為

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

圖1 三針式TDR探頭示意圖

將式（1）、（2）聯立可得［27］：

1.2 數值建模

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

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

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

圖3 含水合物沉積物模型的幾何結構及其網格

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

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

聯立此二方程得：

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

式中，Φ 是標量函數，表示垂直于探針長度方向的截面內的電勢分布，滿足拉普拉斯方程［29］。

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

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

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

圖4 數值模型中三針式TDR探頭與研究區域

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

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

1.3 數據分析方法

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

式中：Ω代表研究區域；Φ0（x，y）是研究區域中具有均勻分布介電常數被測介質時的電勢分布；Φ（x，y）是研究區域中為實際被測介質時的電勢分布。空間權重函數描述了TDR探頭在研究區域Ω內任何點的測量靈敏度。空間權重函數滿足：

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

式中，K（x，y）是研究區域Ω中點（x，y）處被測介質的相對介電常數。

利用研究區域的電勢分布可以計算出探針周圍電場的有效極化面積。有效極化面積是指對TDR 測量響應有主要貢獻區域的面積，在該區域之外，被測介質性質的改變對TDR測量響應沒有顯著的影響［24］。基于空間權重函數的概念，定義垂直探針方向截面上的部分區域對TDR 測量響應的貢獻在總響應中所占的比例f，即

在數值模型中，式（12）中分子為部分有限元網格區域內網格的空間權重函數wi與網格面積Ai乘積的和。從內探針周圍開始，按照空間權重函數從大到小的順序進行排序，whi為研究區域內空間權重函數的最大值。當f＝90%時，認為此時區域面積為TDR的有效極化面積，當f＞90%時，區域面積的增加對TDR測量響應的影響不顯著。

2 TDR探頭參數優化

2.1 探針的個數

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

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

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

對于具有不同探針個數的探頭，以水合物位于區域A為例繪制出探針周圍的電場強度分布，如圖7 所示。由于TDR探頭的內外探針間距相同，由電場強度分布圖可知，隨著探針個數由3 增加至5，高電場覆蓋區域形狀由不規則橢圓形變為以外探針為頂點的三角形和四邊形，有效極化面積逐步增大。

圖7 不同探針個數時水合物位于區域A時電場強度

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

圖8 水合物位于不同區域時測得表觀介電常數的相對變化量與探針個數的關系

2.2 探針的間隔與直徑

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

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

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

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

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

3 水合物飽和度仿真測量

通過改變多孔介質數值模型中水合物區域的大小（見圖3），分別模擬0%～100%水合物飽和度下的含水合物沉積物，并仿真分析水合物飽和度與表觀介電常數的關系。其中，表觀介電常數用式（11）計算得到。

以數值模型的計算結果為基礎，對水合物飽和度Sh與Ka進行擬合得到：

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

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

Topp經驗式：

Wright經驗式：

CRIM模型：

MD模型：

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

圖11 本文仿真擬合與其他經驗式或模型對比圖

4 結語

針對TDR探頭在含水合物沉積物中的應用，基于有限元數值模型計算了垂直于探針方向平面上測試區域內含水合物沉積物的TDR 響應，在對探針個數、探針間距和直徑等結構參數進行優化設計的基礎上，對多孔介質中水合物飽和度進行了仿真測量。得到以下三點認識：

（1）隨著探針個數由3 增加至5，高電場覆蓋區域形狀由不規則橢圓形變為以外探針為頂點的三角形和四邊形，有效極化面積逐步增大；三針式TDR 探頭對于內外探針之間的水合物探測靈敏度最高。

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

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