Fabrication and Application of Drill String Dynamics Experiment Bench for Small-Scale Horizontal Wells
-
摘要:
水平井钻井过程中钻头与钻柱的运动规律非常复杂，通过现场测试研究钻头与钻柱振动的方法成本高，无法准确解耦各因素，且数据采样率低，不能进行系统分析。因此，以物理相似为基础，设计了一套基于微钻头−钻柱−岩石的小尺度水平井钻柱动力学模拟实验台架；在钻柱上设计分段套管，模拟套管对钻柱的约束；为更加真实反应钻头与岩石互作用对钻杆及钻头振动的影响，设计微型PDC钻头进行模拟破岩；开展了水平井钻柱室内试验研究，分析了钻柱系统的耦合振动、转速及钻压对钻柱系统振动的影响规律。结果表明：钻头的轴向振动对钻头的扭转振动有较大的影响，且当钻柱发生扭转振动时，其功率谱密度明显高于匀速转动下钻头的功率谱密度；随着转速增大，钻头和钻柱的横向振动加速度也随之增大，钻头扭转振动在主频上的振动幅值也随之增大；随着钻压增大，钻头粘滑振动和横向振动加剧，但对轴向振动的影响不明显。
Abstract:In the process of horizontal well drilling, the motion law of the drill bit and the drill string is very complex. The cost of studies by conducting site tests for the vibration of the drill bit and the drill string is high, which makes it difficult to accurately decouple various factors. In addition, the data sampling rate is low, hindering the systematic analysis. Therefore, based on physical similarity, a drill string dynamics simulation experiment bench for small-scale horizontal wells was designed on the basis of micro drill bit, drill string, and rock. The segmented casing was designed on the drill string to simulate the constraint of casing on the drill string. In order to more truly reflect the impact of the interaction between drill bit and rock on drill pipe and drill bit vibration, a micro polycrystalline diamond compact (PDC) bit was designed to simulate rock breaking. Laboratory experiments of drill strings in horizontal wells were carried out to study the effect law of coupling vibration, rotational speed, and weight on bit (WOB) on the vibration of the drill string system. The results show that the axial vibration of the drill bit has a great influence on the torsional vibration of the drill bit. When the drill string is under torsional vibration, the power spectral density of the drill bit is obviously higher than that under constant rotation. With the increase in the rotation speed, the transverse vibration acceleration of the drill bit and drill string also increases, and the vibration amplitude of the drill bit under torsional vibration within the main frequency increases. With the increase in WOB, stick-slip vibration and transverse vibration of the drill bit are intensified, but the axial vibration is not significantly affected.
-
钻进裂缝性地层过程中，钻井液漏失是最常见的井下复杂问题之一[1–2]。井漏不仅会使钻井周期增长，处理不当甚至会导致井眼报废[3]。堵漏时一般根据经验判断裂缝开度，优选堵漏材料的配方和粒径[4–5]。研究表明，堵漏材料粒径小于0.6倍缝宽时，堵漏材料可以保持较高浓度进入裂缝内，形成较为致密的封堵层[6]；但实际应用中，漏失层位、漏失速度和钻井液密度基本相同的2口井采用相同的堵漏材料，实际堵漏效果相差很大。裂缝粗糙度对其中流体运移和固体颗粒的铺置有较大影响[7–8]。天然裂缝表面起伏大，在何处形成封堵层、封堵材料粒径如何选配是堵漏需要考虑的核心问题
裂缝性地层的钻井液漏失模型是根据钻井液在裂缝和基质中的流动方程建立的，结合边界条件求解，进而模拟漏失速度与钻井液性质、裂缝形态的关系。F. Sanfillippo等人[9]假设钻井液为牛顿流体，研究了一条无限长裂缝与井筒相交时的钻井液漏失规律。A. Lavrov等人[10–11]建立了有限长单裂缝钻井液漏失模型，定量分析了漏失规律。O. Lietard等人[12]提出了单条光滑裂缝条件下的宾汉模式钻井液漏失模型。R. Majidi等人[13–14]建立了光滑径向裂缝条件下赫巴模式钻井液漏失模型。M. Ozdemirtas等人[15]提出了考虑分形裂缝粗糙度的二维钻井液漏失模型。Xia Yang等人[16–17]建立了二维离散裂缝网络钻井液漏失模型，并采用数值模型研究了裂缝开度与漏失速度的关系。贾利春等人[18–19]引入裂缝力学开度和裂缝迂曲度，表征裂缝面粗糙度对钻井液漏失的影响，分析了承压堵漏效果与堵漏材料的特性。李大奇等人[20–21]基于蒙特卡罗方法生成三维离散裂缝网络，建立了考虑裂缝线性变形的钻井液漏失模型，研究了分形粗糙面对钻井液漏失的影响。吕开河等人[22]建立了径向粗糙裂缝漏失模型；王明波[23]采用连续随机叠加法和中点位移法生成分形裂缝表面，研究了单条粗糙裂缝的漏失规律。截至目前，国内外关于三维裂缝性地层的钻井液漏失模型，尚未考虑粗糙裂缝网络对钻井液漏失速度的影响。因此，笔者在椭圆坐标系下改进了W-M分形函数，建立了三维椭圆裂缝粗糙表面形貌方程，表征了具有分形粗糙特征的椭圆裂缝表面；建立了三维裂缝性地层钻井液漏失流固耦合力学模型，并采用扩展有限单元法对模型进行了数值离散和求解，研究了裂缝表面分形参数对漏失量的影响，揭示了裂缝表面的粗糙度对钻井液漏失的影响规律，为裂缝开度预测及堵漏材料粒径优选提供了理论支撑。
1. 三维椭圆裂缝粗糙表面形貌方程
研究表明，裂缝表面是粗糙的，具有自仿射特征，可以用分形理论研究[21,24]。分形理论可以用来模拟和表征具有分形特征粗糙表面的形貌，常用的方法有随机中点位移法、分形布朗运动法、逐次随机增加法、逆傅里叶变换法和W-M函数法。前4种方法通常会在模拟过程中带入随机变量，并且产生迭代运算，模拟不规则几何形状的表面时有较大的难度；W-M函数是空间和粗糙表面轮廓频率的方程，能够唯一确定已知分形维数和特征尺度参数时的粗糙表面轮廓曲线和形貌。
根据实际观察，大多数天然裂缝面具有椭圆或近似椭圆的形状。Baecher圆盘模型和Veneziano多边形模型都可用于构建裂缝平面模型，但Veneziano模型在三维空间中的几何形状是相当复杂的，因此大部分研究者使用相对简单的 Baecher模型[25]。目前，三维裂缝扩展模型中，椭圆裂缝多采用几何模型[26–27]，因此选择椭圆面来表征裂缝。传统的W-M函数是在直角坐标系下建立起伏面与空间位置的关系，对于生成椭圆粗糙表面具有一定难度，为此引入椭圆坐标系
(ξ,η) ，其与直角坐标系(x,y) 的关系为：x=Lccoshξcosη (1) y=Lcsinhξsinη (2) 式中：Lc为椭圆坐标系下共焦椭圆的焦距，m；
(ξ,η) 是以Lc为焦距的椭圆坐标系下的坐标值，m；(x,y) 是直角坐标系下的坐标值，m。若给定一个椭圆裂缝的长轴和短轴，可求得焦距Lc，进而求得直角坐标系空间坐标所对应的椭圆坐标。M. Ausloos等人[28]给出了直角坐标系下各向同性三维分形表面的W-M函数。在岩石断裂面上，中部的高度高于边缘的高度[29]，因此可以在椭圆坐标系下改进W-M分形函数，并引入特征尺度衰减系数，使粗糙表面轮廓高度由裂缝中央向裂缝边界逐渐减小，改进的W-M函数可以表示为：
z(ξ,η)=(ξe−ξξe)1D−1GD−2(lnγM)12M∑m=1∑n=nmin (3) 式中：z为粗糙表面的高度，m；
{\xi }_{\mathrm{e}} 为椭圆裂缝边界坐标值；D为分形维数，对于三维粗糙表面，2 < D < 3 ；G为分形表面的高度幅值，m；\gamma 为频率密度因子，对于满足正态分布的粗糙表面通常取\gamma = 1.5 [30]；M为构造表面的粗糙峰叠加数量；n为表面上微凸体的频率指数；nmin为所有微凸体中最小的频率指数；nmax为最大频率指数；k0为表面波数，{\mathrm{m}}^{-1} ，与取样长度L之间满足关系k0=2π/L；\varphi m,n为位于（0,2π）内的随机相位。式（3）等号右边第一项为特征尺度衰减系数，反映了椭圆粗糙表面高度由中心向边界衰减，随着
\xi 增大，距裂缝边界距离减小，衰减系数取值减小，裂缝面起伏高度减小。给定椭圆的长轴为100 m，短轴为50 m，表面高度幅值为0.5 mm，微凸体频率指数nmin=5、nmax=50，粗糙峰叠加数量M=20，分形维数D分别取2.2，2.4，2.6和2.8，利用式（3）生成了椭圆裂缝的粗糙表面，结果如图1所示。从图1可以看出，随着分形维数增大，椭圆表面的轮廓变得越密集和复杂，出现了更多的凸峰和凹谷，轮廓幅值也逐渐增大，可见分形维数显著影响椭圆裂缝表面的粗糙程度。
![]()
图 1 改进的W-M函数模拟不同分形维数的椭圆粗糙表面Figure 1. Simulated elliptical rough surfaces with different fractal dimensions by improved WM function裂缝表面的粗糙度主要由分形维数D、高度幅值G和粗糙峰叠加数量M表征。分形维数反映了裂缝表面整体的自相似性程度，高度幅值反映了裂缝表面的起伏程度，粗糙峰叠加数量反映了裂缝表面粗糙凸起的数量，三者共同决定了裂缝表面的粗糙程度。
2. 裂缝性地层钻井液漏失模型的建立
基于多孔弹性力学构建了钻井液漏失流固耦合模型，修正了立方定律，形成了基质与裂缝耦合的钻井液漏失数学模型。
2.1 基质流固耦合方程
基质的应力平衡方程为：
\nabla \cdot{\boldsymbol{ \sigma }}= 0 (4) 柯西总应力张量
{\boldsymbol{\sigma }} 可写为：{\boldsymbol{\sigma}} = {\boldsymbol{\sigma}} ' - \alpha p{\boldsymbol{I}} (5) 式中：
{\boldsymbol{\sigma }}' 为有效应力，Pa（若采用线弹性本构，则有效应力张量与应变张量的关系满足{\boldsymbol{\sigma}} ' = {\boldsymbol{C}}:{\boldsymbol{\varepsilon }} ）；C为线弹性材料的各向同性张量，Pa；应变张量与位移矢量满足关系{\boldsymbol{ \varepsilon}} = \left( {\nabla {\boldsymbol{u}} + {{\left( {\nabla {\boldsymbol{u}}} \right)}^{\rm{T}}}} \right)/2 ；p为孔隙压力，Pa；\alpha 为Biot系数；I为单位张量。以
\delta {\boldsymbol{u}} 为试函数，将式（5）代入式（4），可得变形控制方程的弱形式：\begin{split} &{\int }_{\varOmega }\nabla \delta {\boldsymbol{u}}(C:\varepsilon -\alpha p{\boldsymbol{I}}){\rm{d}}\varOmega -{\int }_{{\varGamma }_{{\rm{q}}}}\llbracket \delta {\boldsymbol{u}} \rrbracket \cdot p{\rm{d}}\varGamma -\\ &\qquad\qquad{\int }_{{\varGamma }_{t}}\delta u\cdot {{\boldsymbol{\tau}}}{\rm{d}}\varGamma =0 \end{split} (6) 式中：
{\varGamma _t} 为应力边界；{\varGamma _q} 为水压边界；{\boldsymbol{\tau}} 为边界应力矢量，Pa。流体在基质中的流动方程为：
\frac{\partial \zeta }{\partial t}+\nabla \cdot {\boldsymbol{q}}=0 (7) 式中：
\zeta 为单位体积的流体体积增量，m3；{\boldsymbol{ q}} 为单位面积的流量，\mathrm{m}/\mathrm{s} 。流体压力与流体体积增量的关系为：
p = Q(\zeta - \alpha e) (8) 式中：
Q 为Biot模量[31]，Pa；e 为体应变。由式（7）和式（8），可得基质流体流动方程为：
\alpha \nabla \cdot \dot{\boldsymbol{ u }}+ \frac{1}{Q}\frac{{\partial p}}{{\partial t}} + \nabla \cdot{\boldsymbol{q}} = 0 (9) 以
\delta p 为试函数，代入达西公式，可得基质流体流动方程的弱形式：\begin{split} {\int }_{{\varGamma }_{{\rm{f}}}}\delta p{q}_{\mathrm{d}}{\rm{d}}\varGamma =&{\int }_{\varOmega }\alpha \delta p\nabla \cdot \boldsymbol{\dot{u}}{\rm{d}}\varOmega +{\int }_{\varOmega }{Q}^{-1}\delta p\dot{p}{\rm{d}}\varOmega +\\&{\int }_{{\varGamma }_{{\rm{q}}}}\delta p \bar q{\rm{d}}\varGamma +{\int }_{\varOmega }\frac{{K}_{\mathrm{m}}}{\mu }\nabla \left(\delta p\right)\nabla p{\rm{d}}\varOmega \end{split} (10) 式中：
{\varGamma _{\rm{f}}} 为裂缝壁面；{{q}}_{\mathrm{d}} 为钻井液由基质向单位面积裂缝的滤失量，\mathrm{m}/\mathrm{s} ；{\varGamma _{\rm{q}}} 为水压边界；\bar q 为定流量，m/s；{{K}}_{\mathrm{m}} 为基质渗透率，{\mathrm{m}}^{2} ；\mu 为流体黏度，Pa·s。2.2 裂缝流动方程
假设裂缝内流体的流动满足立方定律：
{\boldsymbol{q}}_{\mathrm{f}}=-\frac{{w}^{3}}{12\mu }{\nabla }p_{\mathrm{f}} (11) 式中：
{\boldsymbol{q}}_{\mathrm{f}} 为裂缝内单位裂缝长度的截面流量，{\mathrm{m}}^{2}/\mathrm{s} ；w 为裂缝开度，m；pf为裂缝内流体压力，Pa。对于粗糙裂缝，需要对立方定律进行修正。取粗糙表面上任意一段微元，微元在x-y平面的投影长度分别为dx和dy，如果选取的微元足够小，那么微元可以看作由2个光滑平板组成。沿z方向裂缝表观开度为
{w_{\rm{b}}} ，即裂缝上下表面的垂直距离，微元与x与y轴的夹角为ax和ay，微元长度和开度分别为{\rm{d}}\eta 和{\rm{d}}\xi ：{\rm{d}}\eta = \frac{{{\rm{d}}x}}{{\cos {a_x}}} (12) {\rm{d}}\xi = \frac{{{\rm{d}}y}}{{\cos {a_y}}} (13) 设微元两条边的单位向量分别为：
{{\boldsymbol{e}}_\eta } = \left( {\cos {a_x},0,\sin {a_x}} \right) (14) {{\boldsymbol{e}}_\xi } = \left( {0,\cos {a_y},\sin {a_y}} \right) (15) 微元平面的单位法向量
{{\boldsymbol{e}}_n} 为：{{\boldsymbol{e}}_n} = \frac{{{{\boldsymbol{e}}_\eta } \times {{\boldsymbol{e}}_\xi }}}{{\left| {{{\boldsymbol{e}}_\eta } \times {{\boldsymbol{e}}_\xi }} \right|}} (16) 裂缝的实际开度为上下粗糙表面的最短距离，是表观开度在法向量方向上的投影：
w = \left| {{w_{\rm{b}}} \cdot {{\boldsymbol{e}}_n}} \right| = \frac{{{w_{\rm{b}}}\cos {a_x}\cos {a_y}}}{{\sqrt {1 - {{\sin }^2}{a_x}{{\sin }^2}{a_y}} }} (17) 将式（12）、式（13）、式（17）代入式（11），可得粗糙裂缝修正立方定律：
{{\boldsymbol{q}}_{\mathrm{f}}} = - {\boldsymbol{J}}\frac{{w_{\rm{b}}^3}}{{12\mu }}{\nabla}p (18) 其中 {\boldsymbol{J}} = \frac{1}{{{{\left( {1 - {{\sin }^2}{a_x}{{\sin }^2}{a_y}} \right)}^{\frac{3}{2}}}}} \left[ {\begin{array}{*{20}{c}} {{{\cos }^4}{a_x}{{\cos }^2}{a_y}} &0 \\ 0& {{{\cos }^2}{a_x}{{\cos }^4}{a_y}} \end{array}} \right] (19) 裂缝内流体流动的质量守恒方程为：
\frac{{\partial w}}{{\partial t}} + \nabla \cdot {{\boldsymbol{q}}_{\rm{f}}} = {Q_{{\rm{inj}}}}\delta \left( s \right) (20) 式中：
{Q_{\mathrm{inj}}} 为单位截面的裂缝内的源汇项，m/s；\delta \left( s \right) 为狄拉克函数；{w_{\rm{b}}} 为裂缝表观开度，m。以
\delta p 为试函数，可得质量守恒方程的弱形式为：\int\limits_{{\varOmega _{\rm{f}}}} {\delta p\frac{{\partial w}}{{\partial t}}{\rm{d}}\varOmega } + \int\limits_{{\varOmega _{\rm{f}}}} {\delta p\nabla \cdot {{\boldsymbol{q}}_{\rm{f}}}{\rm{d}}\varOmega } - \int\limits_{{\varOmega _f}} {\delta p{Q_{{\rm{inj}}}}\delta (s){\rm{d}}\varOmega } = 0 (21) 其中
{\int }_{{\varOmega }_{{\rm{f}}}}\delta p\nabla \cdot {{\boldsymbol{q}}}_{{\rm{f}}}{\rm{d}}\varOmega ={\int }_{{\varGamma }_{{\rm{f}}}}\delta p{{\boldsymbol{q}}}_{{\rm{f}}}\cdot {\boldsymbol{n}}_{{\rm{f}}}{\rm{d}}\varGamma +{\int }_{{\varOmega }_{{\rm{f}}}}{\nabla }\delta p\boldsymbol{J}\frac{{w}_{{\rm{b}}}^{3}}{12\mu }{\nabla }p{\rm{d}}\varOmega (22) 式中，
{\boldsymbol{n}}_{{\rm{f}}} 为裂缝面的法向量。式（22）等号后第一项为钻井液从裂缝向基质的滤失速度，将式（22）代入式（21），忽略裂缝内的源汇，可得滤失速度的表达式：
{\int }_{{\varGamma }_{{\rm{f}}}}\delta p{{\boldsymbol{q}}}_{\mathrm{f}}\cdot {\boldsymbol{n}}_{{\rm{f}}}{\rm{d}}\varGamma =-{\int }_{{\varOmega }_{{\rm{f}}}}{\nabla }\delta p\boldsymbol{J}\frac{{w}_{{\rm{b}}}^{3}}{12\mu }{\nabla }p{\rm{d}}\varOmega -{\int }_{{\varOmega }_{{\rm{f}}}}\delta p\frac{\partial w}{\partial t}{\rm{d}}\varOmega (23) 3. 扩展有限元数值模拟方法
采用扩展有限元法对钻井液漏失模型进行离散和求解时，通过局部加强形函数捕捉位移和压力特征，将裂缝与基质解耦[32]，可以大幅提升数值计算的精度和效率。不包含裂缝的普通网格节点记为Ns，裂缝完全穿过的网格节点记为Nc，包含裂缝边缘的网格节点记为Nt（见图2，图中蓝色节点为Nc，红色节点为Nt，其他节点为Ns）。
3.1 固体位移加强方式
对于裂缝完全穿过的单元，节点采用Heaviside函数加强：
H(f) = \left\{ {\begin{array}{*{20}{l}} {1,{\text{ }}f \geqslant 0} \\ {0,{\text{ }}f < 0} \end{array}} \right. (24) 式中：f为水平集函数，f >0，表示在裂缝的上侧，f <0表示在裂缝的下侧。
对于缝尖，应力存在奇异性，对于各向同性弹性材料，采用缝尖分支函数加强[32]：
\begin{split} &\qquad\qquad [{F_l}(r,\theta ),{\text{ }}l = 1,2,3,4] = \\ &\left[\sqrt r \sin \frac{\theta }{2}{\text{ }}\sqrt r \cos \frac{\theta }{2}{\text{ }}\sqrt r \sin \frac{\theta }{2}{\text{sin}}\theta {\text{ }}\sqrt r \cos \frac{\theta }{2}\sin \theta \right] \end{split} (25) 式中：r为某点到裂缝边缘的垂直距离，m；θ为x轴正向与r径向的夹角，（°）。
采用平移加强形函数的方法来削弱位移场混合单元，平移修正后的位移逼近式写为：
\begin{split} u(x) = &\sum\limits_{i \in {N_{\rm{s}}}} {{N_i}(x){u_i}} + \sum\limits_{j \in {N_{\rm{c}}}} {{N_j}(x)(H(x) - H({x_j}))} {a_j} +\\& \sum\limits_{k \in {N_{\rm{t}}}} {{N_k}(x)\sum\limits_{l = 1}^4 {({F_l}(x) - {F_l}({x_k}))} } b_k^l = {\boldsymbol{N}} \cdot {\boldsymbol{U}} \end{split} (26) 式中：Ni，Nj和Nk为标准有限元形函数；aj，
b_k^l 为2类加强函数对应的固体位移附加自由度；N为由标准有限元形函数和固体位移加强形函数组成的形函数矢量；U为由标准有限元自由度和位移附加自由度组成的位移自由度矢量。3.2 流体压力加强方式
流体压力在裂缝两侧连续，但流体压力梯度在裂缝两侧不连续。Moës采用改进的绝对值函数加强裂缝完全穿过的网格节点[33]：
{\phi }^{{\rm{c}}}\left(x\right)={\sum }_{j}\left|f\left({x}_{j}\right)\right|{N}_{j}\left(x\right)-\left|{\sum }_{j}f\left({x}_{j}\right){N}_{j}\left(x\right)\right| (27) 式中：
{\phi }^{{\rm{c}}}\left(x\right) 为裂缝完全穿过网格的加强函数；Nj为标准有限元形函数。流体速度在裂缝边缘有奇异性，Xia Yang等人[34]根据裂缝处压力的渐近解析解，构造了加强形函数：
{\phi ^{\rm{t}}}({\boldsymbol{x}}) = \sqrt r \cos \frac{\theta }{2} (28) 式中，
{\phi }^{{\rm{t}}}\left(x\right) 为包含裂缝边缘网格的加强函数。同样，采用平移加强形函数的方法来削弱压力场混合单元，平移修正后的流体压力逼近式写为：
\begin{split} p\left(x\right)=&{\sum }_{i\in {N}_{{\rm{s}}}}{N}_{i}\left(x\right){p}_{i}+{\sum }_{a\in {N}_{{\rm{c}}}}{N}_{a}\left(x\right){\phi }^{{\rm{c}}}\left(x\right){\tilde{p}}_{a}+\\ &{\sum }_{b\in {N}_{{\rm{t}}}}{N}_{b}\left(x\right)({\phi }^{{\rm{t}}}(x)-{\phi }^{{\rm{t}}}({x}_{b})){\tilde{p}}_{b}={\boldsymbol{H}}\cdot {\boldsymbol{P}} \end{split} (29) 式中：Ni，Na和Nb为标准有限元形函数；
{p}_{i} 为标准有限元自由度；{\tilde p_a} 和{\tilde p_b} 为裂缝完全穿过网格与包含裂缝边缘网格加强函数对应的流体压力附加自由度；H为由标准有限元形函数和流体加强形函数组成的压力场形函数矢量；P为由标准有限元自由度和流体压力附加自由度组成的压力自由度矢量。3.3 数值离散方程
根据扩展有限元中位移的不连续加强方式，受力变形后的裂缝开度用不连续自由度计算：
\begin{split} {w}_{\mathrm{f}}=&\llbracket u\rrbracket\cdot \boldsymbol{{n}_{\mathrm{f}}}=\left({\sum }_{i\in {N}_{{\rm{c}}}}{N}_{i}\left(x\right){a}_{i}+2{\sum }_{q\in {N}_{{\rm{t}}}}\sqrt{r}{N}_{q}\left(x\right){b}_{q}\right)\cdot\\ & \boldsymbol{{n}_{\mathrm{f}}}=\llbracket N\rrbracket\cdot \boldsymbol{U}\cdot \boldsymbol{{n}_{\mathrm{f}}} \\[-1pt] \end{split} (30) 将位移和压力逼近表达式代入控制方程积分弱形式，可以得到钻井液漏失模型的离散计算格式：
\left[\begin{array}{cc}0& 0\\ {M}_{\mathrm{u}}& {M}_{\mathrm{p}}\end{array}\right]\left[\begin{array}{c}\dot{U}\\ \dot{P}\end{array}\right]+\left[\begin{array}{cc}{K}_{\mathrm{u}}& {K}_{\mathrm{p}}\\ 0& {K}_{\mathrm{p}\mathrm{p}}\end{array}\right]\left[\begin{array}{c}U\\ P\end{array}\right]=\left[\begin{array}{c}{F}_{1}\\ {F}_{2}\end{array}\right] (31) 其中\quad \left\{\begin{array}{l}{K}_{\mathrm{u}}={\displaystyle \int }_{\varOmega }{\left(\nabla {\boldsymbol{N}}\right)}^{{\rm{T}}}:{\boldsymbol{C}}:\nabla {\boldsymbol{N}}{\rm{d}}\varOmega \\ {K}_{\mathrm{p}}=-{\displaystyle \int }_{\varOmega }{\left(\nabla {\boldsymbol{N}}\right)}^{{\rm{T}}}\alpha {{\boldsymbol{m}}}^{{\rm{T}}}{{\boldsymbol{H}}}^{{\rm{T}}}{\rm{d}}\varOmega +{\displaystyle \int }_{{\varGamma }_{{\rm{f}}}}\llbracket N\rrbracket {{\boldsymbol{H}}}^{{\rm{T}}}{\rm{d}}\varGamma \\ {F}_{1}={\displaystyle \int }_{{\varGamma }_{t}}\boldsymbol{N}\cdot \boldsymbol{\tau}{\rm{d}}\varGamma \end{array}\right. (32) \left\{ \begin{array}{l}{M}_{\mathrm{u}}={\displaystyle \int }_{\varOmega }\alpha {\boldsymbol{Hm}}\nabla {\boldsymbol{H}}{\rm{d}}\varOmega +{\displaystyle \int }_{{\varGamma }_{{\rm{f}}}}H{\llbracket N\rrbracket }^{{\rm{T}}}{\rm{d}}\varGamma \\ {M}_{\mathrm{p}}={\displaystyle \int }_{\varOmega }{Q}^{-1}{\boldsymbol{H}}{{\boldsymbol{H}}}^{{\rm{T}}}{\rm{d}}\varOmega +{\displaystyle \int }_{{\varGamma }_{{\rm{f}}}}\dfrac{w}{{K}_{l}}{\boldsymbol{H}}{{\boldsymbol{H}}}^{{\rm{T}}}{\rm{d}}\varGamma \\ {K}_{\mathrm{p}\mathrm{p}}={\displaystyle \int }_{\varOmega }\dfrac{{k}_{\mathrm{m}}}{\mu }{\left(\nabla {\boldsymbol{H}}\right)}^{{\rm{T}}}\nabla {\boldsymbol{H}}{\rm{d}}\varOmega +{\displaystyle \int }_{{\varGamma }_{{\rm{f}}}}\dfrac{{w}_{\mathrm{b}}^{3}}{12\mu }{\left(\nabla {\boldsymbol{H}}\right)}_{\mathrm{f}}\cdot J\cdot {\left(\nabla {\boldsymbol{H}}\right)}_{\mathrm{f}}^{{\rm{T}}}{\rm{d}}\varGamma \\ {F}_{2}={\displaystyle \int }_{{\varGamma }_{q}}H{\bar q}{\rm{d}}\varGamma \end{array}\right. (33) {\boldsymbol{m}}=[1,0,0] (34) 式中，
{\left(\nabla H\right)}_{\mathrm{f}} 为压力场形函数矢量在裂缝平面的梯度。4. 粗糙裂缝对钻井液漏失量的影响
4.1 物理模型
首先根据粗糙表面形貌方程生成2个椭圆裂缝表面Z1和Z2，假设裂缝初始张开度为w0，此时裂缝内各点处的表观开度wb=Z1–Z2+w0，若wb小于0则设置为0。三维漏失模型如图3所示。模拟地层的长宽高均为200 m，图3中黑线为直井筒，蓝色椭圆面为裂缝，3条裂缝相互沟通形成漏失通道。当井筒内压力高于地层压力时，钻井液由井筒沿着裂缝向地层内部漏失。最小水平主应力沿着X轴，最大水平主应力沿着Y轴。模型基本参数为：地层渗透率5 mD，钻井液黏度10 mPa·s，Biot 系数0.90，Biot 模量10 GPa，井筒压力50 MPa，岩石弹性模量20 GPa，最大水平主应力70 MPa，上覆地层压力75 MPa，裂缝初始张开度2.0 mm，粗糙裂缝高度幅值0.5 mm，粗糙峰叠加数量20，分形维数2.8，地层压力40 MPa，岩石泊松比0.25，最小水平主应力45 MPa。
利用扩展有限元方法求解所提出的裂缝性地层钻井液漏失模型，得到不同漏失时间下的地层压力（见图4），可以看出，随着漏失时间增长，裂缝面的压力明显升高，这是因为裂缝是钻井液漏失的主要通道，而钻井液由裂缝向基质的滤失是一个相对缓慢的过程。图5所示为裂缝上点A（130 m，80 m，105 m）处压力和开度随漏失时间的变化。由图5可知，裂缝的开度在漏失一定时间内不变，这是因为，缝内流体压力小于最小水平主应力，无法将裂缝进一步撑开；随着漏失时间增长，缝内流体压力大于最小水平主应力，裂缝开度开始增大；漏失约50 h后，裂缝开度由初始1.4 mm增大至2.6 mm。
![]()
图 4 不同漏失时间下的地层压力分布Figure 4. Formation pressure distribution with different time of leakage![]()
图 5 裂缝上点A处压力和开度随漏失时间的变化规律Figure 5. Variation of pressure and fracture opening on point A of fracture with leakage time4.2 分形维数对漏失量的影响
分形维数反映了裂缝表面整体的自相似性程度。改变4.1节物理模型的分形维数，其他参数保持不变，得到不同分形维数下钻井液的累计漏失量（见图6）。不同分形维数下裂缝面高斯积分点的最小开度、最大开度、平均开度和接触率（开度为0的高斯积分点所占比例）见表1。
![]()
图 6 不同分形维数下钻井液的漏失量Figure 6. Drilling fluid leakage amount with different fractal dimensions表 1 不同分形维数下裂缝开度的分布Table 1. Distribution of fracture opening with different fractal dimensions分形维数 最小开度/mm 最大开度/mm 平均开度/mm 接触率，% 2.1 1.89 2.12 2.01 0 2.3 1.48 2.66 2.03 0 2.5 0.82 3.27 2.06 0 2.7 0 4.58 2.17 3.4 2.9 0 7.50 2.08 8.2 取不同分形维数，对照表1所列裂缝属性及图6对应的漏失量，分析分形维数对漏失量的影响：随着分形维数增大，平均开度变化不明显，裂缝最小开度减小，裂缝面接触的区域增加（此区域内钻井液不能流动），裂缝最大开度增大，钻井液漏失量明显增大。这说明，分形维度增大，虽然对裂缝的平均开度影响不大，但由于裂缝表面更为粗糙，使裂缝中部分区域的开度显著增大，导致钻井液更容易流动，从而钻井液漏失量更大。
4.3 高度幅值对漏失量的影响
高度幅值反映了裂缝表面起伏程度。改变4.1节所提及物理模型的高度幅值，其他参数保持不变，得到不同高度幅值下钻井液的累计漏失量（见图7）。不同高度幅值下裂缝面高斯积分点的最小开度、最大开度、平均开度和接触率见表2。
![]()
图 7 不同高度幅值下钻井液的漏失量Figure 7. Drilling fluid leakage amount with different height amplitudes表 2 不同高度幅值下裂缝开度分布Table 2. Distribution of fracture opening with different height amplitudes高度幅值/mm 最小开度/mm 最大开度/mm 平均开度/mm 接触率，% 0.1 1.23 2.69 1.98 0 0.2 0.57 3.29 1.96 0 0.5 0 4.94 1.92 0.31 1.0 0 7.50 1.88 5.21 2.0 0 12.25 1.94 16.45 取不同高度幅值，对照表2所列裂缝属性与图7对应的漏失量，分析得到高度幅值加数量对漏失量的影响：随着高度幅值增大，虽然平均开度变化不明显，裂缝的最小开度减小，但最大开度增大，钻井液漏失量明显增大。由于高度幅值反映了裂缝表面的起伏程度，是裂缝表面粗糙度的重要指标之一，所以在裂缝平均开度大致相同的前提下，裂缝表面起伏程度越大，钻井液的漏失量越大。
4.4 粗糙峰叠加数量对漏失量的影响
粗糙峰叠加数量反映了裂缝表面粗糙凸起的数量。改变粗糙峰叠加数量，其他参数保持不变，得到不同粗糙峰叠加数量下钻井液的累计漏失量，结果如图8所示。为进一步分析粗糙峰叠加数量对漏失量的影响，统计了不同粗糙峰叠加数量下裂缝面高斯积分点的最小开度、最大开度、平均开度和接触率，结果见表3。
![]()
图 8 不同粗糙峰叠加数量下钻井液的漏失量Figure 8. Drilling fluid leakage amount with different number of superimposed rough peaks表 3 不同粗糙峰叠加数量下裂缝开度的分布Table 3. Distribution of fracture opening with different number of superimposed rough peaks粗糙峰
叠加数量最小开度/
mm最大开度/
mm平均开度/
mm接触率，
%5 0 9.18 2.12 0.25 10 0 8.42 2.02 1.28 20 0 6.50 1.88 5.21 50 0 4.92 1.71 17.74 100 0 3.56 1.53 32.65 取不同粗糙峰叠加数量，对照表3所列裂缝属性与图8对应的漏失量，分析粗糙峰叠加数量对漏失量的影响：裂缝面的粗糙凸起越多，裂缝的平均开度越小，裂缝面接触率越大，钻井液漏失量越小。当堵漏材料进入裂缝中，此时堵漏材料可看作裂缝面接触的凸起，增加了裂缝面的接触率，降低了裂缝的渗透率，因此能显著降低钻井液的漏失量。
4.5 裂缝大开度区域对漏失量的影响
前述分析发现，钻井液漏失量与裂缝最大开度呈明显的正相关关系。在4.1节参数的基础上，最大开度上限分别设为5，4，3和2 mm计算钻井液漏失量，并与不限制最大开度的情况进行对比，结果如图9所示。从图9可以看出，限制裂缝的最大开度，钻井液漏失量明显降低，说明裂缝的“大开度区域”控制着钻井液在裂缝内的流量，堵漏应重点在开度较大处架桥。
![]()
图 9 不同最大开度限制下钻井液的漏失量Figure 9. Drilling fluid leakage amount with different maximum opening limits5. 实例计算
塔里木油田博孜某井钻至巴什基奇克组出现钻井液漏失，地层已知参数为：地层渗透率23.8 mD，钻井液黏度40 mPa·s，Biot系数0.87，Biot模量23.5 GPa，井筒压力110 MPa，岩石弹性模量40 GPa，最大水平主应力145 MPa，上覆地层压力105 MPa，最小水平主应力100 MPa，地层压力90 MPa，岩石泊松比0.28。根据成像测井可知，该漏失层段有2条裂缝。该井的漏失情况见图10的“实测漏失点”。将已知参数作为模型输入数据，将分形参数和裂缝初始开度作为待调整变量，拟合现场记录的钻井液漏速。若不考虑粗糙裂缝，裂缝初始开度为0.17 mm时，拟合结果如图10中红线所示；若考虑裂缝粗糙度，改变裂缝初始开度和分形参数，当裂缝初始开度为0.14 mm、分形维数为2.5、高度幅值为0.18 mm、粗糙峰叠加数量为25时，拟合结果如图10中蓝线所示。从图10可以看出，考虑裂缝粗糙表面的分形参数后，拟合结果更好，此时裂缝初始开度更小，但漏失速度更大，说明裂缝粗糙面的分形参数对实际漏失速度有较大的影响。应当指出，记录的漏失数据越多，参数的拟合精度越高。
![]()
图 10 模型计算漏速与实测漏速对比Figure 10. Comparison between calculated and measured loss velocity6. 结论与建议
1）在椭圆坐标系下建立三维裂缝粗糙表面形貌方程，修正立方定律，建立了裂缝性地层钻井液漏失流固耦合力学模型，采用扩展有限单元法进行求解，形成了三维裂缝性地层钻井液漏失数学模型。
2）采用W-M分形函数所描述的裂缝表面粗糙度主要由分形维数、高度幅值和粗糙峰叠加数量等3个参数表征，分别反映裂缝表面整体的自相似性程度、起伏程度和裂缝表面粗糙凸起的叠加数量。
3）分形维数和高度幅值不影响裂缝平均开度，但与钻井液漏失量呈正相关关系；粗糙峰叠加数量与裂缝的平均开度和钻井液漏失量呈负相关关系；钻井液漏失量与裂缝最大开度有明显的正相关关系；反演裂缝开度时需要充分考虑裂缝表面粗糙度的影响，堵漏材料应重点在开度较大处架桥。
4）该研究需要进一步考虑钻井液漏失过程对岩石力学参数的影响，评价漏失中裂缝扩展造成漏失量增大的风险。
-
![]()
图 5 微型PDC钻头切削齿分布和刀翼轮廓
Figure 5. Distribution of cutting teeth and profile of the wing on micro PDC bit
![]()
图 6 钻头动力学响应曲线及转速频率分析（钻压=173N）
Figure 6. Dynamics response curve and rotational speed frequency analysis of drill bit （WOB = 173 N）
![]()
图 7 钻头转速频率分析（钻压=173 N）
Figure 7. Rotational speed frequency analysis of drill bit（WOB=173 N）
![]()
图 8 钻头角位移随时间变化及钻头速度相轨迹
Figure 8. Variation of angular displacement of bit with time and speed phase trajectory of drill bit
![]()
图 9 钻头转速及钻柱扭转振动强度级值随时间变化的曲线
Figure 9. Variation of rotation speed of drill bit and torsional vibration level of drill string with time
![]()
图 10 钻头转速与轴向加速度随时间变化的曲线
Figure 10. Variation of rotation speed and axial acceleration of drill bit with time
![]()
图 11 钻头扭转振动与横向振动的关系
Figure 11. Relationship between torsional vibration and axial vibration of drill bit
![]()
图 12 钻头扭转振动功率谱密度与归一化频率的关系
Figure 12. Relationship between power spectral density and normalized frequency of torsional vibration of bit
表 1 微型PDC钻头的布齿参数
Table 1 Parameters of micro PDC bit tooth arrangement
刀翼号 齿号 径向位置/mm 轴向高度/mm 方位角/（°） 法向角/（°） 侧倾角/（°） 后倾角/（°） 齿长度/mm 1 1 4.00 34.87 0 －20.00 －3 －15 5 2 1 7.26 35.17 120 －20.00 －3 15 5 3 1 10.51 35.47 240 －20.00 －2 15 5 1 2 13.77 35.87 0 －17.35 －1 10 5 2 2 16.94 36.18 120 5.77 －1 10 5 3 2 20.07 35.18 240 29.71 －1 10 5 1 3 22.66 32.84 0 52.53 －1 5 5 2 3 24.44 29.66 120 68.95 －1 5 5 3 3 25.44 26.21 240 73.83 －1 5 5 1 4 26.00 22.82 0 90.00 0 0 5 
下载: 导出CSV
表 2 转速与钻头扭转振动的关系
Table 2 Relationship between rotation speed and torsional vibration of drill bit
转速/
（r·min−1）钻压/
N驱动频率/
Hz主频/
Hz基频/
Hz转速最大振幅/
（r·min−1）30 177 0.50 12.2 0.50 4.0 45 145 0.75 11.8 0.73 5.8 60 176 1.00 11.8 1.00 6.8 75 144 1.25 12.32 1.23 7.3
下载: 导出CSV
表 3 钻压与钻头轴向、横向振动的关系
Table 3 Relationship between WOB，and axial and transverse vibration of drill bit
钻压/
N扭转振动
强度轴向加速度/（m·s−2） 横向加速度/（m·s−2） 均值 极大值 均值 极大值 149 1.17 0.036 0.080 0.350 2.060 180 1.35 0.036 0.075 0.640 2.930 212 1.43 0.036 0.077 0.810 3.160 242 1.57 0.036 0.080 0.890 3.680
下载: 导出CSV
-
[1] 王敏生,光新军,耿黎东. 页岩油高效开发钻井完井关键技术及发展方向[J]. 石油钻探技术,2019,47(5):1–10. WANG Minsheng, GUANG Xinjun, GENG Lidong. Key drilling/completion technologies and development trends in the efficient development of shale oil[J]. Petroleum Drilling Techniques, 2019, 47(5): 1–10.
[2] QIU Hongyuan, YANG Jianming, BUTT S, et al. Investigation on random vibration of a drillstring[J]. Journal of Sound and Vibration, 2017, 406: 74–88. doi: 10.1016/j.jsv.2017.06.016
[3] 高德利,黄文君. 超深井工程理论与技术若干研究进展及发展建议[J].石油钻探技术,2024,52(2):1-11. GAO Deli, HUANG Wenjun. Research and development suggestions on theory and techniques in ultra-deep well engineering[J]. Petroleum Drilling Techniques, 2024, 52(2): 1-11.
[4] ELSAYED M A, PHUNG C C. Modeling of drillstrings[C]//ASME 2005 24th International Conference on Offshore Mechanics and Arctic Engineering. New York: ASME, 2005: 15-24.
[5] NAVARRO-LÓPEZ E M. An alternative characterization of bit-sticking phenomena in a multi-degree-of-freedom controlled drillstring[J]. Nonlinear Analysis: Real World Applications, 2009, 10(5): 3162–3174. doi: 10.1016/j.nonrwa.2008.10.025
[6] LIAO C M, BALACHANDRAN B, KARKOUB M, et al. Drill-string dynamics: reduced-order models and experimental studies[J]. Journal of Vibration and Acoustics, 2011, 133(4): 041008. doi: 10.1115/1.4003406
[7] 姚永汉,狄勤丰,朱卫平,等. 底部钻具组合的涡动特征分析[J]. 石油钻探技术,2010,38(4):84–88. YAO Yonghan, DI Qinfeng, ZHU Weiping, et al. Analysis of whirling properties of bottom hole assembly[J]. Petroleum Drilling Techniques, 2010, 38(4): 84–88.
[8] 石祥超,焦烨,刘景涛,等. 考虑深井井下动力钻具影响的钻柱粘滑振动规律[J]. 天然气工业,2024,44(6):87–97. SHI Xiangchao, JIAO Ye, LIU Jingtao, et al. Stick-slip vibration law of drill strings under the influence of mud motors in deep wells[J]. Natural Gas Industry, 2024, 44(6): 87–97.
[9] 吕苗荣,沈诗刚. 钻柱黏滑振动动力学研究[J]. 西南石油大学学报(自然科学版),2014,36(6):150–159. LYU Miaorong, SHEN Shigang. The simulation and analysis of drillstring stick-slip vibration[J]. Journal of Southwest Petroleum University(Science & Technology Edition), 2014, 36(6): 150–159.
[10] 滕学清,狄勤丰,李宁,等. 超深井钻柱粘滑振动特征的测量与分析[J]. 石油钻探技术,2017,45(2):32–39. TENG Xueqing, DI Qinfeng, LI Ning, et al. Measurement and analysis of stick-slip characteristics of drill string in ultra-deep wells[J]. Petroleum Drilling Techniques, 2017, 45(2): 32–39.
[11] WILSON J K, WHITACRE T, HEISIG G. High-frequency at-bit measurements provide new insights into torsional dynamics when drilling with steerable mud motors in unconventional horizontal wells[R]. SPE 194072, 2019.
[12] CHEN Shilin, WISINGER J, DUNBAR B, et al. Identification and mitigation of friction- and cutting action-induced stick-slip vibrations with PDC bits[R]. SPE 199639, 2020.
[13] RAYMOND D W, ELSAYED M A, POLSKY Y, et al. Laboratory simulation of drill bit dynamics using a model-based servohydraulic controller[J]. Journal of Energy Resources Technology, 2008, 130(4): 043013.
[14] KHULIEF Y A, AL-SULAIMAN F A. Laboratory investigation of drillstring vibrations[J]. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2009, 223(10): 2249–2262. doi: 10.1243/09544062JMES1550
[15] LU Haochuan, DUMON J, CANUDAS-DE-WIT C. Experimental study of the D-OSKIL mechanism for controlling the stick-slip oscillations in a drilling laboratory testbed[C]//2009 IEEE Control Applications, (CCA) & Intelligent Control, (ISIC). Piscataway, NJ: IEEE Press, 2009: 1551-1556.
[16] FRANCA L F P. Drilling action of roller-cone bits: modeling and experimental validation[J]. Journal of Energy Resources Technology, 2010, 132(4): 043101. doi: 10.1115/1.4003168
[17] FORSTER I, MACFARLANE A, DINNIE R. Asymmetric vibration damping tool-small scale rig testing and full scale field testing[R]. SPE 128458, 2010.
[18] FORSTER I. Axial excitation as a means of stick slip mitigation-small scale rig testing and full scale field testing[R]. SPE 139830, 2011.
[19] VLAJIC N, LIAO C M, KARKI H, et al. Draft: stick-slip motions of a rotor-stator system[J]. Journal of Vibration and Acoustics, 2014, 136(2): 021005. doi: 10.1115/1.4025994
[20] 张东霄,李波,邱琳琳,等. 基于相似原理的钻柱振动模拟试验研究[J]. 矿山机械,2016,44(5):8–12. ZHANG Dongxiao, LI Bo, QIU Linlin, et al. Simulation test study on vibration of drill string based on similarity principle[J]. Mining & Processing Equipment, 2016, 44(5): 8–12.
[21] 温欣,管志川,邵冬冬,等. 大斜度井眼中钻柱运动特性模拟试验[J]. 中国石油大学学报(自然科学版),2018,42(2):79–86. doi: 10.3969/j.issn.1673-5005.2018.02.009 WEN Xin, GUAN Zhichuan, SHAO Dongdong, et al. Experimental study of dynamic characteristics of drillstring in highly deviated well[J]. Journal of China University of Petroleum (Edition of Natural Science), 2018, 42(2): 79–86. doi: 10.3969/j.issn.1673-5005.2018.02.009
[22] 任福深,刘建华,李洋,等. 基于相似理论的水平受压钻柱试验装置研制[J]. 石油机械,2016,44(4):12–16. REN Fushen, LIU Jianhua, LI Yang, et al. Development of similarity theory-based horizontal well drill string test apparatus[J]. China Petroleum Machinery, 2016, 44(4): 12–16.
[23] LIAN Zhanghua, ZHANG Qiang, LIN Tiejun, et al. Experimental and numerical study of drill string dynamics in gas drilling of horizontal wells[J]. Journal of Natural Gas Science and Engineering, 2015, 27(part 3): 1412-1420.
[24] WESTERMANN H, GORELIK I, RUDAT J, et al. A new test rig for experimental studies of drillstring vibrations[J]. SPE Drilling & Completion, 2015, 30(2): 119–128.
[25] LI Wei, HUANG Genlu, NI Hongjian, et al. Modeling and experimental study on motion states of laboratory-scale bottom hole assembly in horizontal wells[J]. Energies, 2020, 13(4): 925. doi: 10.3390/en13040925
[26] ABDO J, HASSAN E, AL-SHABIBI A. Fabrication of experimental setup for investigation of drill string fatigue failure[C]//ASME 2017 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. New York: ASME, 2017: V008T12A055.
[27] ESMAEILI A, ELAHIFAR B, FRUHWIRTH R K, et al. Laboratory scale control of drilling parameters to enhance rate of penetration and reduce drill string vibration[R]. SPE 160872, 2012.
[28] 邹德永,王瑞和. 刀翼式PDC钻头的侧向力平衡设计[J]. 中国石油大学学报(自然科学版),2005,29(2):42–44. doi: 10.3321/j.issn:1000-5870.2005.02.010 ZOU Deyong, WANG Ruihe. Lateral force balancing design of blade PDC bits[J]. Journal of China University of Petroleum(Edition of Natural Science), 2005, 29(2): 42–44. doi: 10.3321/j.issn:1000-5870.2005.02.010
[29] XUE Qilong. Data analytics for drilling engineering: theory, algorithms, experiments, software[M]. Cham: Springer, 2020.
-
期刊类型引用(3)
1. 杨万有，郑春峰，李昂，王磊. 海上电泵结蜡井热循环洗井工艺参数优化设计. 西南石油大学学报(自然科学版). 2019(01): 129-136 . 
百度学术
2. 朱广海，刘章聪，熊旭东，宋洵成，王军恒，翁博. 电加热稠油热采井筒温度场数值计算方法. 石油钻探技术. 2019(05): 110-115 . 本站查看
3. 沈园园. 南堡潜山高温油气藏井下循环温度的数值模拟. 断块油气田. 2017(04): 570-573 . 百度学术
其他类型引用(2)
计量
- 文章访问数: 234
- HTML全文浏览量: 50
- PDF下载量: 78
- 被引次数: 5
