基于有限元模型的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实验测试，可以获得不同水合物饱和度条件下的表观介电常数实测数据，进而综合利用有限元数值模拟和实验测试数据来建立可靠性高、准确度高、适用范围广的水合物饱和度计算模型。