石油与天然气化工  2024, Vol. 53 Issue (1): 77-82
微芯片技术在油气田井下漏层定位中的应用
苗海龙1,2 , 马跃2 , 石钊睿3 , 邱正松1     
1. 中国石油大学(华东)石油工程学院;
2. 中海油田服务股份有限公司油田化学事业部;
3. 浙江探芯科技有限公司
摘要目的 准确地定位油气田开发过程中的井漏位置,为制定堵漏措施提供支持,提高堵漏的成功率。方法 由于漏失层位大量钻井液漏失会造成该位置温度变化异常,通过一种搭载微芯片的微型井下测量球状仪在井内循环,获取井下温度、压力分布数据,然后基于这些数据掌握井下温度分布规律,并与理想的温度分布规律进行对比,从而比较准确地定位漏失点位置。结果 通过捕获井口返出的微芯片,获取井内温度压力数据,然后再对温度梯度变化数据进行平滑处理,最终能够基于数据分析出较为准确的漏失层位。结论 基于微芯片测量数据来定位漏失点是一种针对井下漏失情况诊断的新型、低成本且便捷的方法,为创新型井下测量微芯片的应用提供了一个非常有价值和代表性的场景。
关键词微芯片    温度分布    漏点判断    井漏    
Application of microchip technology in the location of downhole circulation loss layer in oil and gas fields
MIAO Hailong1,2 , MA Yue2 , SHI Zhaorui3 , QIU Zhengsong1     
1. School of Petroleum Engineering, China University of Petroleum (East China), Qingdao, Shandong, China;
2. Department of Oilfield Chemistry, China Oilfield Services Limited., Yanjiao, Hebei, China;
3. Zhejiang Tanxin Technology Co., Ltd., Ningbo, Zhejiang, China
Abstract: Objective The purpose of this study is to accurately locate the circulation loss layer during the drilling of oil and gas fields, and provide support for the control of loss, improve the success rate of curing loss. Methods Due to the abnormal temperature change caused by a large amount of drilling fluid loss in the circulation loss layer, a downhole measuring ball equipped with a microchip was circulated in the well to obtain downhole temperature and pressure distribution data, and then the downhole temperature distribution law was established based on these data, compared with the ideal temperature distribution law, which can more accurately locate the circulation loss layer. Results The experiment shows that the temperature and pressure data in the well can be obtained by capturing the microchips returned from the wellhead, and then the temperature gradient change data can be smoothed, and the more accurate loss zone can be analyzed by using the data. Conclusions Locating leakage points based on microchip measurement data is a new, low-cost and convenient method for downhole leakage diagnosis, which is a very valuable and representative application scenario of innovative downhole measurement microchips.
Key words: microchip    temperature distribution    location of the circulation loss layer    lost circulation    

在油气田开发过程中,井漏是一种后果较为严重的复杂情况,可导致钻井液大量损失,如果井漏失控,会发生井控和井壁垮塌的严重情况,从而造成非生产时间增加和巨大的经济损失[1]。对于堵漏作业来说,首先需要解决的问题是漏失点的定位。过去的文献中提到了几种定位漏失地层的方法,但这些方法需要先另外下入井下测量工具以获取目标地层的信息,然后再定位漏失点,从而大大增加作业的成本。所以,研究一种经济性好的准确定位漏失层位的方法是非常重要的。

针对上述问题,研究人员开展了一系列研究。Shea[2]开展了漏失点深度定位的研究,通过在钻井液中加入放射性的材料后将其泵入井下,然后下入录井仪器追踪放射信号来定位漏失钻井液的位置。Owings[3]提到了一种通过对比实际循环压力场和正常时循环压力场来定位漏失地层的方法。Han介绍了一种使用多普勒效应的原位测量泥浆流量装置。Dacoord等[5]介绍了一种管状的装置,其外部安装有传感器,用于测量流体的流速并将数据传输至地表,而这个数据可以用于定位漏失层和分析漏失的严重程度。Dakin等[6]阐述了一种使用4种不同液体来定位漏失层的方法。Horne[7]在现场通过注入冷却水测试来确定漏失点,在该方法中,地表的冷却水被泵入井下,同时下入录井工具进行温度采集,观察温度变化以判断漏点。Dokhani等[8]开发了一种瞬态热交换模拟算法来计算井筒循环时和静止时的压力场和温度场情况。该算法考虑了包含泥浆漏失量及排量变化等多种现场参数的影响。

井下测量微芯片是一种微尺寸的球形或胶囊形仪器,内部包含微控制单元、内存、传感器、数据传输模块及微型可充电锂电池,所有的电路和元器件被包裹在球形或胶囊形的保护性材料内,整个系统的直径约为8~12 mm。微芯片所使用的保护性材料具有极高的强度和玻璃化温度,可以保证其内部电路及元器件能承受井下的高温高压环境。微芯片可以在井下随钻井液循环并回收,测量并存储流经位置的温度和压力数据,回收后可以下载读取这些数据并进行分析和处理。

朱祖杨研究了一种基于微芯片示踪器实现漏失位置检测的方法:首先下入多个测量短节,然后再投入一批微芯片示踪器进行测量,测量完之后再通过几个不同位置的测量短节数据对微芯片数据进行深度标定,最后通过使用标定后的数据观察温度突变点,从而确定漏失点。

1 理论模型

基于研究使用微芯片测量数据来定位漏失点及确定漏速,其主要原理是通过钻井液发生漏失后引起钻井液流速的改变,进而导致整个温度场温度发生变化,通过微芯片所测量数据来监控该变化,进而定位其发生变化所在的层位。主要的理论方法为通过建立井下的瞬态传热模型来进行分析。以某深水井为例,其井筒的结构和钻井液的流动情况如图 1所示(1 ℉=17.2 2℃,1 ft=0.305 m)。在漏失发生之前,钻井液以一个固定的流量(wp)流经钻杆内和环空段。假设在环空段有一个轴向的流体均匀通过,那么钻井液的流速会因为分流至有漏失的地层而降低。而由于钻井液流速降低,整个井下钻杆、环空和地层之间的热交换平衡将发生改变。同时,由于流量降低,流经环空漏失点上端的钻井液将会被进一步冷却。该技术既可以应用在陆地井,也可以应用在海上井。但应用于海上井时,要考虑海水对温度场的影响。由于海上井较陆地井分析更为复杂,在模型推导中,为了更清晰地进行理论分析,故以海上井为例,把整个井下结构分为3段:第1段为海床以上部分;第2段为海床到漏失点之间的部分;第3段为漏失点以下部分。

图 1     与环境温度梯度对应的井下结构分段图

井下的热量可以源于一团物质的进入而带来的焓,也可来自外界物质的热传导、热对流及热辐射。同样的,热量也可通过相同的方式离开系统。累积的整个能量都是以内部能量形式存在的(井下的动势能可忽略)。系统整体的热平衡见式(1)~式(4):

$ \frac{\mathrm{d} \varPhi}{\mathrm{d} t}=\sum\nolimits_{k=1}^M \dot{m}_k \hat{H}_k+\sum\nolimits_{j=1}^N \dot{Q}_j $ (1)

式中:Φ为单位体积的内能,BTU(1 BTU= 1.056 kW);$\dot{m}_k$为质量流量,lb/h(1 lb/h=0.45 kg/h);$\hat{H}_k$为比焓,Btu/lb;$\dot{Q}_j$为热传递速率,Btu/h。

$ \mathrm{d} \varPhi=\left(\frac{\partial \varPhi}{\partial t}\right)_V \mathrm{~d} T+\left(\frac{\partial \varPhi}{\partial V}\right)_T \mathrm{~d} V $ (2)

式中:V为体积,ft3(1 ft3=0.028 m3);T为时间, h。

如果不考虑体积因为时间推移发生的变化,那么式(2)中的第2项可以被省略,即:

$ \mathrm{d} \varPhi=\left(\frac{\partial \varPhi}{\partial t}\right)_V \mathrm{~d} T=C_V m \mathrm{~d} T $ (3)

式中: CV为单位体积的比热,BTU/(lb·℉);m为质量,lb;T为时间,h。

结合式(1)和式(3),可得式(4):

$ \hat{H}=\frac{\varPhi+P V}{m}=u+\frac{P}{\rho} $ (4)

式中:P为压力,lb/ft2(1 lb/ft2=47.88 Pa);u为单位内能,Btu/lb;ρ为密度,lb/ft3(1 lb/ft3=16 kg/m3)。

在通常的井下热传导应用中,式(4)中的第1项比内能u远远大于式中第2项的值。因此,式(4)中的第2项在实际应用中可以省略。对井筒的3段分别单独应用热平衡进行分析,从而可以分别获取钻井液在钻杆内、环空及地层的支配方程。具体推导出的热平衡方程见式(5):

钻杆内(第1段、第2段、第3段):

$ \frac{\rho_{\mathrm{m}} c_{\mathrm{m}} A_{\mathrm{p}}}{2 {\rm{\mathsf{π}}} r_{\mathrm{p}} U_{\mathrm{ap}}} \frac{\partial T_{\mathrm{p}}}{\partial \mathrm{t}}=-\frac{\dot{m}_{\mathrm{p}} c_{\mathrm{m}}}{2 {\rm{\mathsf{π}}} r_{\mathrm{p}} U_{\mathrm{ap}}} \frac{\partial T_{\mathrm{p}}}{\partial z}+\left(T_{\mathrm{a}}-T_{\mathrm{p}}\right) $ (5)

式中:ρm为泥浆密度,lb/ft3cm为泥浆比热,Btu/(lb·℉);Ap为钻杆横截面积,ft2(1 ft2=0.093 m2);rp为钻杆半径,ft(1 ft=0.3 m);Uap为钻杆与环控之间的传热系数,Btu/(hr·℉· ft2);Tp为钻杆流体温度,℉(℉=32+℃×1.8);$\dot{m}_{\mathrm{p}}$为钻杆质量流量,lb/hr;Ta为环空流体温度,℉;z为井的垂直深度,ft。

环空(第1段):

$ \begin{aligned} \frac{\rho_{\mathrm{m}} c_{\mathrm{m}} A_{\mathrm{a}}}{2 {\rm{\mathsf{π}}} r_{\mathrm{w}}^R h_{\mathrm{ar}}^R} \frac{\partial T_{\mathrm{a}}^R(z, t)}{\partial_{\mathrm{t}}}= & \frac{c_{\mathrm{m}} \dot{m}_{\mathrm{a}}^U}{2 {\rm{\mathsf{π}}} r_{\mathrm{w}}^R h_{\mathrm{ar}}^R} \frac{\partial T_{\mathrm{a}}^R(z, t)}{\partial z}- \\ & \left(T_{\mathrm{a}}^R-T_{\mathrm{w}}\right)+\frac{r_{\mathrm{P}} U_{\mathrm{ap}}}{r_{\mathrm{w}}^R h_{\mathrm{ar}}^R}\left(T_{\mathrm{a}}-T_{\mathrm{a}}^R\right) \end{aligned} $ (6)

式中: Aa为环空横截面积,ft2rwR为立管段的井眼半径,ft;harR为立管和环空之间的对流传热系数,Btu/(hr·℉·ft2);TaR为立管段的环空流体温度,℉;$\dot{m}_{\mathrm{a}}^U$为漏失点上段的环空质量流量,lb/hr;Tw为井筒壁温度,℉。

环空(第2段):

$ \begin{aligned} \frac{\rho_{\mathrm{m}} c_{\mathrm{m}} A_{\mathrm{a}}}{2 {\rm{\mathsf{π}}} r_{\mathrm{w}} h_{\mathrm{af}}} \frac{\partial T_{\mathrm{a}}^u(z, t)}{\partial_{\mathrm{t}}}= & \frac{c_{\mathrm{m}} \dot{m}_{\mathrm{a}}^U}{2 {\rm{\mathsf{π}}} r_{\mathrm{w}} h_{\mathrm{af}}} \frac{\partial T_{\mathrm{a}}^u(z, t)}{\partial \mathrm{z}}-\left(T_{\mathrm{a}}^u-T_{\mathrm{w}}\right)+ \\ & \frac{r_{\mathrm{P}} U_{\mathrm{ap}}}{r_{\mathrm{w}} h_{\mathrm{af}}}\left(T_{\mathrm{p}}-T_{\mathrm{a}}^u\right) \end{aligned} $ (7)

式中: rw为立管段的井眼半径,ft;haf为环空流体和地层之间的对流传热系数, Btu/(hr·℉·ft2);Tau为漏失点上段的环空流体温度,℉。

环空(第3段):

$ \begin{aligned} \frac{\rho_{\mathrm{m}} c_{\mathrm{m}} A_{\mathrm{a}}}{2 {\rm{\mathsf{π}}} r_w h_{\mathrm{af}}} \frac{\partial T_{\mathrm{a}}^u(z, t)}{\partial_{\mathrm{t}}}= & \frac{c_{\mathrm{m}} \dot{m}_{\mathrm{a}}^U}{2 {\rm{\mathsf{π}}} r_{\mathrm{w}} h_{\mathrm{af}}} \frac{\partial T_{\mathrm{a}}^u(z, t)}{\partial z}-\left(T_{\mathrm{a}}^u-T_{\mathrm{w}}\right)+ \\ & \frac{r_{\mathrm{P}} U_{\mathrm{ap}}}{r_{\mathrm{w}} h_{\mathrm{af}}}\left(T_{\mathrm{p}}-T_{\mathrm{a}}^u\right) \end{aligned} $ (8)

式中: haf为环空流体和地层之间的对流传热系数,Btu/(hr·℉·ft2)。

立管(第1段):

$ \dot{Q}_{\mathrm{ar}}=2 {\rm{\mathsf{π}}} r_{\mathrm{r} 0}^R \Delta z h_{\mathrm{r} 0}^R\left(T_{\mathrm{ro}}-T_{\mathrm{a}}^R\right) $ (9)

式中: $\dot{Q}_{\mathrm{ar}}$为立管与环空流体之间热传递速率,Btu/hr;rr0R为立管外侧半径,ft;hr0R为立管外侧热传递系数,Btu/(hr·℉·ft2);Tro为立管外表面温度,℉;TaR为立管段环空流体温度,℉。

式(9)被定义为立管的边界条件外部表面温度,因此,在计算中可被假设为正常的海水温度分布:

$ T_{\mathrm{ro}}=T_{\mathrm{i}}^R(z)=\left\{\begin{array}{l} 5 \cdot 10^{-6} z^2-0.0201 z+59.945, 0 \leqslant z \leqslant 3000 \\ 40, 3000<z \leqslant 5000 \end{array}\right. $ (10)

地层(第2段、第3段):

环空段的泥浆和地层之间的热交换通过井壁发生,因此可得:

$ 2 {\rm{\mathsf{π}}} r_{\mathrm{w}} h_{\mathrm{af}}\left(T_{\mathrm{f}}\left(z, r=r_{\mathrm{w}}, t\right)-T_{\mathrm{a}}\right)=2 {\rm{\mathsf{π}}} r_{\mathrm{w}} k_{\mathrm{f}} \frac{\partial T_{\mathrm{f}}\left(z, r=r_{\mathrm{w}}, t\right)}{\partial r} $ (11)

式中:Tf为地层流体温度,℉;kf为地层的热导率,Btu/(h·ft·℉)。

地层的热交换主要来自热传导和热对流,支配方程见式(12):

$ \begin{aligned} \frac{\partial T_{\mathrm{f}}(z, r, t)}{\partial t}= & \alpha_{\mathrm{f}}\left(\frac{\partial^2 T_{\mathrm{f}}(z, r, t)}{\partial r^2}+\frac{1}{r} \frac{\partial T_{\mathrm{f}}(z, r, t)}{\partial r}\right)+ \\ & \frac{\rho_{\mathrm{fl}} c_{\mathrm{fl}}}{\rho_{\mathrm{f}} c_{\mathrm{f}}} \frac{\partial P_{\mathrm{f}}(z, r, t)}{\partial r} \frac{\partial T_{\mathrm{f}}(z, r, t)}{\partial r} \end{aligned} $ (12)

式中:αf为地层热扩散率,ft2/hr;ρf为地层密度,lbm/ft3ρfl为空隙流体密度,lbm/ft3cf为地层比热,Btu/(lbm·℉);cfl为孔隙流体比热,Btu/(lbm·℉);Pf为地层流体压力,lbf/ft2

井下的孔隙压力分布为:

$ \frac{\partial P_{\mathrm{f}}(z, r, t)}{\partial r}=\alpha_{\mathrm{c}}\left[\begin{array}{l} \left(\frac{\partial^2 P_{\mathrm{f}}(z, r, t)}{\partial r^2}+\frac{1}{r} \frac{\partial P_{\mathrm{f}}(z, r, t)}{\partial r}\right)+ \\ C_{\mathrm{fl}}\left(\frac{\partial P_{\mathrm{f}}(z, r, t)}{\partial r}\right)^2 \end{array}\right] $ (13)
$ \alpha_{\mathrm{c}}=\frac{K_{\mathrm{f}}}{\mu_{\mathrm{fl}} \varphi^2 c_{\mathrm{t}}} $ (14)

式中:αc为压力扩散常数,ft2/hr;Kf为地层的热导率,Btu/(h·ft·℉);φ为内能,BTU;μfl为孔隙流体单位内能,BTU。

2 计算方法

通过微芯片现场测试,获取全井筒温度和压力原始测量值,使用深度标定计算方法将原始数据转换为井深数据,并利用移动平均法进行数据平滑处理。在实际计算中,可以通过调整移动平均的时期个数来实现合适的数据平滑处理结果。最后,计算温度梯度(dT/dZ)变化对应井深的完整结果,形成温度梯度(dT/dZ)曲线用于实际的漏点定位。

在实际的漏点判断中,计算获取的环空温度梯度曲线(dT/dZ)一般遵循以下规则:当dT/dZ < 0时,钻井液向上运动时温度逐渐升高(通常发生在环空下段);当dT/dZ>0时,钻井液向上运动时温度逐渐降低(通常发生在环空上段)。当温度梯度(dT/dZ)曲线发生明显的正向突变时,说明了环空中的液体在加速冷却,或是在减速升温,表明可能存在以下情况:①井下结构发生明显变化;②地温梯度发生明显变化;③漏失情况发生。在实际应用中,井下结构的变化和地温梯度的变化情况可以被排除,因此可以确定漏失点位置。

3 基于微芯片测量数据定位漏失点的循环管路实验
3.1 循环管路实验装置

为了进一步验证上述方法,在实验室搭建了一套专门用于微芯片的模拟井筒循环管路实验系统,同时在管路处设置一处开路,系统可以控制此开路来模拟井下漏失的发生。该系统可以很好地观察和记录微芯片在钻杆内和环空的运动状态。整个循环管路长11 m,斜度为75°,系统原理图如图 2所示。

图 2     微芯片循环管路系统原理图

在整个深度标定实验管路系统的基础上,本实验在井深7 m处井筒上额外添加了一处裂缝,整体由哈夫节(Hough section)包围,并接上排水装置,井底有加热装置,最高可加热至50 ℃,最大压力为0.14 MPa。管路系统内安装有多个扶正器,以确保钻杆无偏心问题。

哈夫节又称柔性伸缩器,是一种管道堵漏器。主要结合塑形接触密封和弹性接触密封两种技术,能够在管道表面紧固并贴合,起到密封作用,以控制介质的泄漏。

3.2 实验步骤

(1) 由于循环管路长度有限,测试前需调整微芯片采样频率至最高频率(0.1 Hz)。

(2) 提前打开电子加热装置,并打开泵开启液体循环。

(3) 实验前将微芯片充满电后激活,并开始计时。

(4) 先关闭漏失点水阀,投入激活的微芯片循环测试。

(5) 将模拟漏失点水阀调至半开状态,投入激活的微芯片循环测试。

(6) 将模拟漏失点水阀调至全开状态,投入激活的微芯片循环测试。

(7) 分别记录微芯片在钻杆内和环空段的运动时间。

(8) 微芯片排出至泥浆罐中,回收微芯片,下载数据。

在设计微芯片之初,考虑了井下温度和压力对其传感系统的损害,故采用的材料都是耐高温和高压力的。在实际应用中存在着回收率的问题,只要在振动筛中收集到微芯片并采集其中数据,就意味着应用成功,对于失返性漏失,通过投掷微芯片的方法是无法完成的,另一种方式是利用井下工具携带的微芯片释放信号的方式进行,并通过泥浆脉冲返回数据,与井下测量相比,此方式可以大幅节省成本。

3.3 实验结果和分析

测试完成后,对回收的微芯片进行数据下载,其原始温度测量结果如图 3所示。

图 3     微芯片原始测量数据

分别在3次不同漏失量(无漏失、轻微漏失、漏失)的微芯片原始测量数据的基础上进行深度标定计算,得到3种不同漏失情况下的微芯片温度场数据,结果如图 4所示。

图 4     温度场数据斜率(温度梯度)变化对比

将无漏失、轻微漏失及漏失3种场景下的微芯片温度场数据进行平滑处理后进行对比,可以发现温度场的斜率(温度梯度)在井深7 m附近有明显不同,即无漏失情况下无明显斜率突变,轻微漏失情况下有一定斜率突变,漏失情况下有明显斜率突变,如图 5所示。

图 5     经过深度标定的两枚微芯片测量的温度场数据对比

这个温度梯度的突变说明了在井深7 m处,流体的温度场被加速冷却,在排除了井身结构和外界温度的影响后可以得出结论:井深7 m处附近有漏失点,这与实际的实验装置情况完全一致。实验结果很好地证明了基于微芯片测量数据定位漏失点的理论计算方法的准确性。在更多的现场实际应用案例中,该方法可以准确地提供现场实际漏失点位置。

4 应用测试

测试井是位于沙特一口井深为10 843 ft(3 356 m)直井,该井发生漏失时,进行微芯片测试,目的是判断漏失位置。通过对获取到的两组微芯片(分别命名为微芯片1、微芯片2)的温度数据进行分析,深度标定计算后得到温度场测量数据如图 5所示。两组微芯片温度场数据具有较好的一致性,相互差别在1 ℃以内,结果表明其测量的准确性,同时,发现钻井液最高温度发生在井深大约9 400 ft(2 865 m)的位置,而不是井底。基于以上两组温度测量数据,计算得到温度梯度原始和平滑后的结果,分别如图 6图 7所示。

图 6     微芯片1的温度梯度数据对比

图 7     微芯片2的温度梯度数据对比

图 6图 7显示的结果可以看出,经过平滑处理后的两组微芯片温度变化梯度数据的噪声更小,经过处理后的温度梯度数据平均的波动幅度是0.005 ℃/ft。更重要的一点是,两组微芯片的数据均显示出在裸眼段的3 109~3 200 m(10 200~10 500 ft)处温度梯度有巨大变化,说明在此处有很大可能发生了漏失,这与钻井日报中显示的3 207 m(10 522 ft)上方遇到漏失的结论非常一致。另外,在实际微芯片现场测试完成后,根据现场记录,现场工程师对该井进行了堵漏作业,重点目标堵漏区域为井深3 100~3 200 m的裸眼井段,经过有针对性的堵漏作业,最终顺利实现堵漏并继续钻进。这说明,通过微芯片测量数据预测的漏失点位置确实与实际情况高度匹配,具有很好的准确性。

4 结论

(1) 通过建立井下的瞬态传热模型进行分析,建立了微芯片定位漏点的理论模型,通过分析微芯片的测量结果可以定位漏失点。同时,通过使用压力微芯片的测量数据还可以计算出实际的漏速,这个方法特别需要全井压力场在未发生漏失前的基准值。

(2) 完成基于研究的微芯片深度标定方法的实验室测试,通过对比深度标定计算出的结果和实际测量结果,证明了该深度标定方法的精确性和可靠性,也验证了基于微芯片测量数据进行深度标定和定位漏失点技术的精度和可靠性。

参考文献
[1]
陈乐亮. 国外井漏处理技术发展趋势综述[J]. 石油与天然气化工, 1994, 23(4): 231-237.
[2]
SHEA J L. Method of locating a zone of lost circulation: US2749444A[P]. 1956-06-05.
[3]
OWINGS A J. Method for locating the depth of a drill string washout or lost circulation zone: US4346594[P]. 1982-08-31.
[4]
OWINGS A J. Method for locating the depth of a drill string washout or lost circulation zone: US4346594[P]. 1982-08-31.
[5]
DACCORD G, GUILLOT D. Identification of lost circulation zones: US13/044633[P]. 2011-03-10.
[6]
DAKIN E, MI LLC. Methods for minimizing fluid loss to and determining the locations of lost circulation zones: US13/128456[P]. 2011-09-15.
[7]
HORNE R N. Characterization, evaluation, and interpretation of well data[M]//DIPIPPO R. Geothermal Power Generation. Duxford: Woodhead Publishing, 2016: 141-163.
[8]
DOKHANI V, MA Y, YU M J. Determination of equivalent circulating density of drilling fluids in deepwater drilling[J]. Journal of Natural Gas Science and Engineering, 2016, 34: 1096-1105.