基于COMBIN40单元等效带库仑阻尼的长冲抽油杆柱下行滞后仿真

王岗; 岳欠杯; 许萍文; 陈重; 刘巨保; 赵志伟; 聂曦暕

doi:10.12404/j.issn.1671-1815.2506898

科学技术与工程 ›› 2026, Vol. 26 ›› Issue (23) :9912-9927 DOI: 10.12404/j.issn.1671-1815.2506898

石油、天然气工业

基于COMBIN40单元等效带库仑阻尼的长冲抽油杆柱下行滞后仿真

王岗 ¹ ,
岳欠杯 ¹^,^* ,
许萍文 ² ,
陈重 ³ ,
刘巨保 ¹ ,
赵志伟 ¹ ,
聂曦暕 ¹

作者信息 +

Simulation on Downward Hysteresis of Long-stroke Sucker Rod Based on COMBIN40 Element Equivalent Band with Coulomb Damping

Author information +

文章历史 +

摘要

为了使油田稳产,针对长冲程抽油杆在采油过程中易发生滞后,从而降低采油效率的现象,引入质量-滑块模型,通过不同时刻的解析解与COMBIN40单元等效计算的数值解对比,两者误差率均在0.1%以下,表明COMBIN40单元能准确地描述带库仑阻尼的长冲程抽油杆运动问题。考虑抽油杆在三维空间中的几何非线性及与油管接触非线性特征,引入COMBIN40单元来离散长冲程抽油杆柱;利用COMBIN40单元集多物理效应的优势,可在单一单元内同步表征抽油杆柱的轴向惯性、弹性变形与摩擦特性;并且仅需定义轴向平动自由度,大幅降低整体方程维度及计算复杂度,解决了常规用BEAM188单元计算效率低的问题。通过Newmark法对时域进行离散,采用节点载荷和节点位移循环迭代求解方法,以实际生产参数为边界条件,通过计算得到了地面示功图,并与实测示功图进行了比较,两者十分接近。对大庆油田8口大斜度长冲程低冲次井中的抽油杆柱运行参数进行仿真计算,得到了不同冲程、冲次、采出液黏度、泵径、井斜角等参数作用下抽油杆的运行参数,结果表明:随着冲程、冲次、采出液黏度、泵径、井斜角的增加,抽油杆下行的摩阻力数值增大,发生滞后现象越明显。此结果可为现场长冲程抽油杆下行滞后分析提供理论依据。

Abstract

In order to maintain stable oilfield production, aimming the easy occurrence of the long-stroke sucker rod hysteresis during oil extraction, which reduced production efficiency, a mass-slider model was introduced. Analytical solutions were compared with numerical solutions by using the equivalent COMBIN40 element at different time points. The error was less than 0.1% between the two results, indicating that the COMBIN40 element accurately described the motion of the long-stroke sucker rod with Coulomb damping. Considering the geometric nonlinearity of the sucker rod string in three-dimensional space and the contact nonlinearity between the rod and the tubing, the COMBIN40 element was introduced to discretize the long-stroke sucker rod string. The axial inertia, elastic deformation, and frictional characteristics of the sucker rod string were represented simultaneously within a single element by taking advantage of the multi-physical coupling capability for the COMBIN40 element. Only the axial freedom of motion needed to be defined, which significantly reduced the dimensionality of the global equations and the computational complexity. This approach resolved the low computational efficiency problem by using the conventional BEAM188 element. The time domain was discretized using the Newmark method. A cyclic iteration method based on nodal loads and displacements was employed to solve the system. Using actual production parameters as boundary conditions, the surface dynamometer card was obtained through simulation and compared with the measured dynamometer card. The calculated and measured results were found to be in close agreement. The operating parameters of sucker rod strings was simulated in eight high-deviation, long-stroke, low-frequency wells of the Daqing Oilfield. The motion parameters of the sucker rods were obtained under varying stroke lengths, stroke frequencies, produced fluid viscosities, pump diameters and well deviation angles. The results show that, as the stroke, stroke frequency, produced fluid viscosity, pump diameter, and well deviation angle increased, the downstroke friction of the sucker rods also increased, and the hysteresis effect became more pronounced. These findings provide a theoretical basis for analyzing the downward hysteresis of long-stroke sucker rods in the field.

Graphical abstract

关键词

抽油杆滞后 / 长冲程 / 低冲次 / 下行速度 / 累计摩阻力 / 大斜度井

Key words

sucker rod hysteresis / long-stroke / low-frequency / downward velocity / cumulative friction resistance / high slope well

引用本文

引用格式 ▾

王岗, 岳欠杯, 许萍文, 陈重, 刘巨保, 赵志伟, 聂曦暕. 基于COMBIN40单元等效带库仑阻尼的长冲抽油杆柱下行滞后仿真[J]. 科学技术与工程, 2026, 26(23): 9912-9927 DOI:10.12404/j.issn.1671-1815.2506898

WANG Gang, YUE Qian-bei, XU Ping-wen, CHEN Zhong, LIU Ju-bao, ZHAO Zhi-wei, NIE Xi-jian. Simulation on Downward Hysteresis of Long-stroke Sucker Rod Based on COMBIN40 Element Equivalent Band with Coulomb Damping[J]. Science Technology and Engineering, 2026, 26(23): 9912-9927 DOI:10.12404/j.issn.1671-1815.2506898

登录浏览全文

4963

注册一个新账户忘记密码

随着石油开采的不断进行,中国油田已普遍进入高含水、低产量的开发中后期,稳产增产面临严峻挑战^[1]。在此背景下,长冲程抽油系统因其低频、大行程的运行特征,成为提升排液量、降低杆管交变应力、延长杆柱寿命、实现节能降耗的关键技术方向,在中国各大油田(如大庆、胜利等)得到推广应用。然而,长冲程采油需要的泵筒长度更长,导致泵段的摩阻力增大^[2],尤其是在高含水期杆管润滑恶化的情况下,抽油杆下行过程中易出现井口与井底的速度不一致和位移时滞现象,这一滞后效应不仅导致系统效率下降、能耗上升,还可能引发抽油杆振动带来额外的动载荷,严重影响系统稳定性,甚至造成抽油杆卡滞、检泵周期缩短和生产成本攀升。

针对上述问题,中外的研究已从多个角度展开,形成了涵盖新材料、新方法、新模型与新工艺的全方位技术积累,为提升有杆抽油系统的可靠性、效率与寿命积累了理论与实践基础。在抽油杆新材料应用方面,为同时应对高矿化度采出液导致的电化学腐蚀与交变载荷下的疲劳强度问题,碳纤维复合材料抽油杆被成功提出,并进行了探索性应用,这类材料凭借其优异的耐腐蚀性、高比强度与抗疲劳特性,为在苛刻腐蚀环境中延长杆柱寿命提供了新材料^[3];在抽油杆寿命预测与安全评估方面,智能算法与精细力学分析相结合的趋势日益明显,一方面,赵岩龙等^[4]创新性地将长短时记忆网络(long short-term memory,LSTM)这一深度学习模型引入抽油杆剩余寿命预测,通过挖掘历史工况数据中的时序规律,实现了较传统模型更高的预测精度,为预防性维护提供了数据驱动的新工具;另一方面,针对具体的失效模式,研究不断走向微观与定量化,康红兵等^[5]聚焦于塔河油田常见的抽油杆接箍偏磨问题,通过细致的力学性能测试与承载能力分析,评估了不同磨损深度对螺纹连接可靠性的影响,得到严重磨损是诱发脱扣事故的关键因素的结论;姚飞等^[6]则从作业安全角度出发,针对大杆径杆管起下过程中的防喷隐患,通过结构优化设计与数值模拟,建立了一套完整的防喷脱接器安全性评价方法,提升了现场作业的风险管控能力。在失效机理分析方面,李杰等^[7]对新疆油田火驱区块D级抽油杆的腐蚀失效进行了全面的金相组织与理化性能分析,得到了抽油杆失效为CO₂/H₂S协同腐蚀与高矿化度采出液引起的垢下腐蚀共同作用的结论,为预防抽油杆腐蚀提供了理论支持。周宇等^[8]则针对稠油开采的特殊工况,从流体力学基础理论出发,建立并求解管内液体动力学方程,提出了抽油杆柱黏滞摩擦载荷的计算方法,深化了对高黏度流体环境中杆柱下行阻力来源的认识。在动力学仿真与振动分析方面,建模的维度与精度在不断提升,研究已从早期的纵向振动分析,发展到综合考虑纵向、横向、扭转振动的耦合动力学模型^[9],并进一步集成了井筒温度场变化对杆柱力学特性的影响^[10],建立了更为完善的热固耦合振动模型。在求解方法上,有限元法、有限差分法以及四阶龙格-库塔法等数值技术被广泛应用,显著增强了对杆柱复杂动态行为(包括滞后现象)的模拟与示功图预测能力^[11],这些工作为系统理解抽油杆柱在井下的真实运动状态提供了强大的理论工具。在工艺优化与现场实践方面,研究成果正快速转化为提升生产效益的直接生产力,付强等^[12]开发的稠油井开采优化软件,能够基于模型计算结果智能推荐抽汲参数,实现了从经验决策到科学优化的转变;孔令路^[13]提出“上快下慢”变冲次抽汲技术,通过巧妙调节上下行程速度来减缓偏磨与滞后,在数十口井的成功应用证明了其工程实用价值;Yang等^[14]将抽油杆动态特性研究与蒸汽辅助重力泄油(steam assisted gravity drainage,SAGD) 这一特定开采工艺相结合,其模型预测结果与现场数据高度吻合,为特殊工艺下的有杆泵系统设计提供了理论指导。此外,一些前沿性研究也在不断拓展该领域的认知边界,张砚雯等^[15]通过数值模拟与试验验证相结合,深入揭示了复合材料抽油杆的损伤演化机制与失效形式,为工程应用中的安全性评价提供了理论支撑。Marino等^[16]对干摩擦受迫振动阻尼系统的影响进行了理论研究,为理解抽油系统中库伦摩擦的非线性动力学效应提供了理论指导。

上述研究针对常规井抽油杆滞后与失效防治方面成果显著,然而对特殊的大斜度、长冲程、低冲次这一特殊且日益增多的井况,抽油杆柱滞后问题的系统研究尚且少见。此类井况中杆管接触非线性更强、摩擦效应更为复杂,传统建模与分析方法面临极大挑战。

为此,现聚焦于大斜度长冲程低冲次抽油杆柱下行滞后问题,摒弃传统有限元模型中复杂的面-面或点-面接触定义与反复搜索过程,采用COMBIN40单元来模拟抽油杆杆柱在井眼中的运动状态,该单元内置库伦摩擦模型,能够直接模拟抽油杆-油管及柱塞-泵筒间的非线性滑动摩擦,避免了接触算法带来的高昂计算成本与收敛性问题。利用COMB-IN40单元集轴向惯性、弹性变形与内置摩擦于一体的多物理场特性,将复杂的杆管相互作用简化为单元层面的本构行为,可极大地降低整体动力学方程的维度与非线性程度,使得模型更专注于滞后现象的核心物理机制。内置库伦摩擦模型可直接参与全局等效节点力向量的组装,保证动力学求解过程的数值稳定性与精度,结合Newmark时间积分法进行时域离散,构建一套高效、稳健的数值分析模型,能够快速完成长冲程、低冲次过程的大规模时程计算。基于所建模型,不仅能够精确量化井口悬点与井底柱塞的位移滞后量及速度差,更进一步系统性地探讨泵径、冲程、冲次、流体黏度、井斜角等关键工程参数对抽油杆运动状态及滞后效应的影响规律,可为这类特殊井的抽汲参数优化与杆柱设计提供理论指导。

1 算例验证

建立一个具有库仑阻尼的弹簧-质量系统模型,如图1所示,滑块的质量M=4.54 kg,滑块向右移动了距离Δ(0.025 m)并释放,滑动摩阻力F=μMgcosθ,其中μ为摩阻系数,取μ=0.19,g为重力加速度,θ滑块与滑动表面间的角度,(°)。F在滑块与表面之间起限制滑动的作用,方向随运动方向的改变而改变,计算不同时刻t的位移u,并把计算结果与理论结果做对比。

如图1(a)所示为弹簧-质量系统的物理模型,现用图1(b)所示的有限元模型进行模拟计算,COMBIN40单元为一种组合单元,将滑块与弹簧并联连接。滑块弹簧劲度系数K₁可任意选定,要足够大来减少弹性接触效应,同时还要保证实际可行的积分时间步长,本例积分时间步长为0.202 5/405=0.000 5 s,K₁=1 751 267.72 N/m。弹簧劲度系数K₂=5 253.80 N/m,滑动摩阻力F=8.34 N。计算结果如图2和表1所示。

由图2和表1可知,在t=0.09 s时刻,有限元模型计算得到的位移为0.871 5 m,理论求得的位移为0.872 1 m,两者误差率为0.07%;在t=0.102 s时刻,有限元模型计算得到的位移为0.831 9 m,理论求得的位移为0.831 3 m,两者误差率为0.07%;在t=0.183 s时刻,有限元模型计算得到的位移为-0.748 4 m,理论求得的位移为-0.748 7 m,两者误差率为0.04%;理论求得的滑动摩阻力为8.34 N,有限元计算结果为8.34 N,误差率为0。综上可知,用COMBIN40单元等效计算得到不同时刻的位移结果与理论结果几乎一样,表明COMBIN40单元能准确地描述带库伦阻尼的长冲程抽油杆运动问题。

2 长冲程抽油杆仿真力学模型

抽油杆在大斜度井中主要以轴向运动为主,其结构如图3所示,将抽油杆看作一维轴向连续体,在建立油管内长冲程抽油杆力学模型前做如下假设:①只考虑抽油杆沿其轴线方向的运动,忽略扭转和弯曲;②普通Y系列异步电机特性较稳定,载荷变化对电机转速影响不明显,因此忽略载荷变化导致的电机转速波动^[17];③每级抽油杆具有恒定的横截面积A、密度ρ和弹性模量E;④杆的横向尺寸远小于纵向尺寸,且只考虑轴向位移;⑤研究对象为定向井,并且处于高含水阶段,长冲程抽油杆运行速度缓慢,因此忽略液体摩阻力的影响。

选取井底到井口长冲程抽油杆为研究对象,并考虑油管的约束作用,建立了如图4所示的动力学分析模型。

(1)载荷:抽油杆在下行过程中主要受自重W_R(t)、浮力F_FL(t)、惯性载荷I_NE(t)、扶正器与油管间的摩阻力F_FY(t)、杆管摩阻力F_RT(t)、柱塞与泵筒间的摩阻力F_ZB等影响。

(2)初始条件:抽油杆柱初始状态为自然静止,并处于下死点位置,即t=0时刻的初始条件为

(1)

$t={t}_{0}=0,\delta \left(t\right)=0,\stackrel{·}{\delta }\left(t\right)=0,\stackrel{··}{\delta }\left(t\right)=0$

式(1)中:$\delta \left(t\right)$为t₀时刻的位移。

(3)边界条件:油管固定,抽油杆上端与钢丝绳连接并随着钢丝绳的移动而上下移动。

(4)井口边界:施加一个正弦位移S(t)为

(2)

$S\left(t\right)=\frac{L\left(t\right)}{2}\left[\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{2\omega }{60/N}-\frac{\pi }{2}\right)+1\right]$

式(2)中:L(t)为抽油杆举升的距离,m;ω为抽油杆振动圆频率,Hz;N为抽油杆冲次,次/min。

(5)井底边界:抽油泵柱塞与泵筒间产生滑动摩阻力F_ZB,抽油杆上行时滑动摩阻力方向向下,抽油杆下行时方向向上。井底柱塞为已知力边界条件F_P(t),其数值可根据泵结构及运行参数确定^[18],即

(3)

${F}_{\mathrm{p}}\left(t\right)=\left\{\begin{array}{l}{F}_{\mathrm{u}}\left(t\right)\\ {F}_{\mathrm{d}}\left(t\right)\end{array}\right.$

式(3)中:F_u(t)、F_d(t)分别为上、下冲程中柱塞所受的阻力,其计算公式分别为

(4)$ F_{\mathrm{u}}(t)=\rho_{\mathrm{L}} g L_{\mathrm{P}}(t)\left(A_{\mathrm{P}}-A_{\mathrm{r}}\right)-\rho_{\mathrm{L}} g L_{\mathrm{m}}(t) A_{\mathrm{P}}+F_{\mathrm{ZB}}$

(5)$ F_{\mathrm{d}}(t)=-\rho_{\mathrm{L}} g L_{\mathrm{P}}(t) A_{\mathrm{r}}-F_{\mathrm{g} \mathrm{lzl}}(t)-F_{\mathrm{ZB}}$

式中:ρ_L为液体密度,kg/m³;A_P为柱塞的横截面积,m²;A_r为与泵相连的抽油杆柱横截面积,m²;g为重力系数,N/kg;L_P(t)为泵挂深度,m;L_m(t)为沉没度,m;F_ZB(t)为柱塞与泵筒间的滑动摩阻力,N,${F}_{\mathrm{Z}\mathrm{B}}\left(t\right)=0.94\frac{{D}_{P}}{\delta \left(t\right)}-140$,其中D_P为柱塞直径,m,$\delta \left(t\right)$为柱衬间隙,m;F_glzl(t)为过流阻力,N,其大小为F_glzl(t)=$ \frac{1.5 n_{\mathrm{k}}}{728 C^{2}} A_{\mathrm{P}}^{3} \frac{1+A_{0} / A_{\mathrm{P}}}{A_{0}^{2}} \frac{S(t)^{2} N^{2}}{\rho_{\mathrm{L}} g}$,其中n_k为游动阀个数,C为流量系数,A₀为游动阀孔的断面积,m²,S(t)为冲程,m,N为冲次,次/min。

3 长冲程抽油杆有限元及动力学模型

3.1 长冲程抽油杆有限元模型

建立长冲程抽油杆有限元模型,采用COMBIN40单元对抽油杆进行网格划分,如图5所示。将复杂的连续杆柱离散为一系列质量-弹簧-滑块系统,每个单元代表一根杆柱。抽油杆、扶正器、接箍与油管及柱塞与泵筒间的接触用COMBIN40单元内置库伦摩擦模型来描述。

COMBIN40单元可在单一单元内同时实现抽油杆柱轴向弹性变形(弹簧特性)、运动惯性(质量特性)与非线性摩擦(滑块特性)的模拟,无需额外建立单元间耦合约束,避免拆分建模中“质量-刚度-摩擦”数据传递的累积误差,保证了轴向响应精度。同时COMBIN40单元较需考虑多自由度的单元(如BEAM188),仅需定义沿井深方向(Z方向)的单一自由度,单元的方程维度从6h×6h降低到h×h(h为节点数),从而保证了计算效率。

3.2 COMBIN40单元内置库伦摩擦模型

将每一根抽油杆等效为一个COMBIN40单元,其滑块系统即为内置库伦摩擦模型,COMBIN40单元结构^[19],如图6所示。

COMBIN40单元的内置库仑摩擦模型分为“静态库仑摩擦模型”与“动态库仑摩擦模型”两部分,模拟抽油杆从静止到运动的摩擦状态,这样单元的内置库仑模型就为两折线模型,其摩擦力与位移的关系如图7所示。

由图7可知,当间隙GAP≠0时,代表间隙未闭合,即抽油杆与油管间未产生接触,这时抽油杆与油管间的摩擦力为0;当间隙GAP=0时,即抽油杆与油管产生接触,又分为两种情况:静止与滑动。①当外部激励载荷小于最大静摩擦力时,抽油杆与油管不滑动;②当外部激励载荷大于最大静摩擦力时,抽油杆与油管滑动。

COMBIN40单元的内置库仑摩擦模型为

(6)$ \left\{\begin{array}{l}\delta_{\mathrm{E}}(t)=u_{j}(t)-u_{i}(t)+\mathrm{GAP} \\F_{\text {fric }}(t)=\left\{\begin{array}{c}F_{1}(t)+F_{2}(t)=\left(K_{1}+K_{2}\right) \delta_{\mathrm{E}}(t) \\F_{1}(t)+F_{2}(t) \leqslant F_{\mathrm{s}} \\\mu_{\text {cube }} F_{\text {cube }} \operatorname{sgn}\left[\dot{\delta}_{\mathrm{E}}(t)\right], \\F_{1}(t)+F_{2}(t)>F_{\mathrm{s}}\end{array}\right.\end{array}\right.$

式(6)中:F₁(t)为弹簧1的弹力,N;F₂(t)为弹簧2的弹力,N;F_s为最大静摩擦力,N;u_j(t)为t时刻j节点的位移,m;u_i(t)为t时刻i节点的位移,m;${\delta }_{\mathrm{E}}\left(t\right)$为t时刻两节点的相对位移,m;${\stackrel{·}{\delta }}_{\mathrm{E}}\left(t\right)$为t时刻两节点的相对速度,m/s;sgn为符号函数;K₁为弹簧1的弹簧刚度,N/m;K₂为弹簧2的弹簧刚度,N/m;μ_cube抽油杆与油管间的摩擦系数;F_cube抽油杆与油管间的接触力,N,F_cube=M_rodgsin(θ_i),其中θ_i为任意位置处的井斜角,(°),M_rod为单根抽油杆的质量,kg;g为重力常数,N/kg。

3.3 抽油杆仿真动力学模型

为得到抽油杆每次运动时抽油杆的悬点载荷,故在计算时将单根抽油杆的质量等效在COMBIN40单元在j节点,局部坐标系下抽油杆单元质量矩阵$\boldsymbol{M}^{\mathrm{E}}$为

(7)$ \boldsymbol{M}^{\mathrm{E}}=\rho g L_{\mathrm{dg}} \frac{d^{2}}{4}\left[\begin{array}{ll}0 & 0 \\0 & 1\end{array}\right]$

式(7)中:ρ抽油杆的密度,kg/m³;d为抽油杆杆径,m;L_dg为单根抽油杆长,m。

局部坐标系下抽油杆单元的非线性刚度矩阵$\boldsymbol{K}_{0}^{\mathrm{E}}$为

(8)$\boldsymbol{K}_{0}^{\mathrm{E}}=k\left[\begin{array}{cc}1 & -1 \\-1 & 1\end{array}\right]$

式(8)中:k为抽油杆刚度系数,当抽油杆未发生滑动时k=K₁+K₂,抽油杆发生滑动时k=K₂,K₁在保证弹簧与滑块不发生弹性接触效应,可以任意取值;K₂=EA/L_dg,其中E为抽油杆的弹性模量,Pa,A为抽油杆的截面积,m²。

根据达朗贝尔原理^[20],引入惯性力就可将动力问题转化为相应的静力问题,再用最小势能原理推导动力学方程,得到局部坐标系下抽油杆单元仿真动力学模型为

(9)$\boldsymbol{M}^{\mathrm{E}} \ddot{\delta}_{t}^{\mathrm{E}}+\boldsymbol{K}_{0}^{\mathrm{E}} \delta_{t}^{\mathrm{E}}+\boldsymbol{F}_{t}^{\mathrm{E}, \text { fric }}=\boldsymbol{F}_{t}^{\mathrm{E}, \text { ext }}$

式(9)中: $\boldsymbol{M}^{\mathrm{E}}$、$\boldsymbol{K}_{0}^{\mathrm{E}} $、$\boldsymbol{F}_{t}^{\mathrm{E}, \text { fric }}$、$\boldsymbol{F}_{t}^{\mathrm{E}, \text { ext }}$分别为局部坐标系下的单元质量矩阵、非线性刚度矩阵、摩阻力列阵、外力等效节点力列阵;$\ddot{\delta}_{t}^{\mathrm{E}}$、$\delta_{t}^{\mathrm{E}}$分别为抽油杆单元t时刻的加速度及位移。

对局部坐标系下抽油杆单元的单元质量矩阵、非线性刚度矩阵、摩阻力列阵、外力等效节点力列阵,按节点编号进行叠加,即

(10)$\left\{\begin{array}{ll}\boldsymbol{M}=\sum_{n=1}^{h} \boldsymbol{M}^{\mathrm{E}}, & \boldsymbol{K}=\sum_{n=1}^{h} \boldsymbol{K}_{0}^{\mathrm{E}} \\\boldsymbol{F}_{\text {fric }}=\sum_{n=1}^{h} \boldsymbol{F}_{\text {fric }}^{\mathrm{E}}, & \boldsymbol{F}_{\text {ext }}=\sum_{n=1}^{h} \boldsymbol{F}_{\text {ext }}^{\mathrm{E}}\end{array}\right.$

式(10)中:M、K、$\boldsymbol{F}_{\text {fric }}$、$\boldsymbol{F}_{\text {ext}}$分别为整体的质量矩阵、非线性刚度矩阵、摩阻力列阵、所有外载荷形成的等效节点力列阵;h为抽油杆单元节点数。

经坐标变换及单元拼装得到抽油杆在整体坐标系下的仿真动力学模型为

(11)$\boldsymbol{M} \ddot{\delta}_{t}+\boldsymbol{K} \delta_{t}+\boldsymbol{F}_{t}^{\text {fric }}=\boldsymbol{F}_{t}^{\text {ext }}$

式(11)中: $\ddot{\delta}_{t}$、${\delta}_{t}$分别为抽油杆系统在t时刻的加速度及位移。

4 模型求解

式(11)为关于抽油杆单元位移$\delta \left(t\right)$的变系数二阶偏微分方程组,无法求得方程的解析解,为此,用有限元法将抽油杆离散为h个单元,h+1个节点,节点编号由井口的编号1逐渐增大到井底的编号h+1,采用Newmark-β积分法^[21]与有限元法相结合的方法对式(11)进行离散求解,基于Newmark-β积分法对时间t进行离散,得

(12)$\left\{\begin{array}{l}\dot{\delta}_{t+\Delta t}=\dot{\delta}_{t}+\left[(1-\alpha) \ddot{\delta}_{t}+\alpha \ddot{\delta}_{t+\Delta t}\right] \Delta t \\\delta_{t+\Delta t}=\delta_{t}+\dot{\delta}_{t} \Delta t+\left[\left(\frac{1}{2}-\beta\right) \ddot{\delta}_{t}+\beta \ddot{\delta}_{t+\Delta t}\right] \Delta t^{2}\end{array}\right.$

式(12)中:α、β为根据积分的精度和稳定性可以调整的参数,当$\alpha \ge 0.5$,$\alpha \ge {(\beta +0.5)}^{2}$积分格式为无条件稳定,本文中取α=0.5,β=0.25。

由式(12)得到${\stackrel{·}{\delta }}_{t+\mathrm{\Delta }t}、{\stackrel{··}{\delta }}_{t+\mathrm{\Delta }t}$的表达式为

(13)$\left\{\begin{array}{l}\dot{\delta}_{t+\Delta t}=a_{1}\left(\delta_{t+\Delta t}-\delta_{t}\right)+a_{4} \dot{\delta}_{t}+a_{5} \ddot{\delta}_{t} \\\ddot{\delta}_{t+\Delta t}=a_{0}\left(\delta_{t+\Delta t}-\delta_{t}\right)-a_{2} \dot{\delta}_{t}-a_{3} \ddot{\delta}_{t} \\a_{0}=\frac{1}{\beta \Delta t^{2}}, a_{1}=\frac{\alpha}{\beta \Delta t}, a_{2}=\frac{1}{\beta t} \\a_{3}=\frac{1}{2 \beta}-1, a_{4}=1-\frac{\alpha}{\beta}, a_{5}=\left(1-\frac{\alpha}{2 \beta}\right) \Delta t\end{array}\right.$

在t+Δt时刻,动力学方程[式(11)]变为

(14)$\boldsymbol{M} \ddot{\boldsymbol{\delta}}_{t+\Delta t}+\boldsymbol{K} \delta_{t+\Delta t}+\boldsymbol{F}_{t+\Delta t}^{\text {fric }}=\boldsymbol{F}_{t+\Delta t}^{\text {ext }}$

将式(12)和式(13)代入式(14)中可得

(15)$\tilde{\boldsymbol{K}} \delta_{t+\Delta t}=\tilde{\boldsymbol{F}}_{t+\Delta t}$

式(15)中: $\tilde{\boldsymbol{K}} $、$\tilde{\boldsymbol{F}}$分别为有效刚度矩阵和有效载荷矩阵,其表达式为

(16)$\left\{\begin{array}{l}\tilde{\boldsymbol{K}}=\boldsymbol{K}+a_{0} \boldsymbol{M} \\\tilde{\boldsymbol{F}}_{t+\Delta t}=\boldsymbol{F}_{t+\Delta t}+\boldsymbol{M}\left(a_{0} \delta_{t}+a_{2} \dot{\delta}_{t}+a_{3} \ddot{\delta}_{t}\right)\end{array}\right.$

求解式(15)时通过残差来判断结果是否稳定收敛,如果收敛则输出t+Δt时刻的位移${\delta }_{t+\mathrm{\Delta }t}$,再将求得的${\delta }_{t+\mathrm{\Delta }t}$代入式(13)就可求得抽油杆柱在整体坐标系下,t+Δt时刻运动的速度${\stackrel{·}{\delta }}_{t}$,加速度${\stackrel{··}{\delta }}_{t}$,直至计算到$t+\mathrm{\Delta }t\ge 5T$(计算5个周期),最后输出井口位移${\delta }_{\mathrm{t}\mathrm{o}\mathrm{p}}$,井底位移${\delta }_{\mathrm{b}\mathrm{o}\mathrm{t}\mathrm{t}\mathrm{o}\mathrm{m}}$,井口速度${\stackrel{·}{\delta }}_{\mathrm{t}\mathrm{o}\mathrm{p}}$,井底速度${\stackrel{·}{\delta }}_{\mathrm{b}\mathrm{o}\mathrm{t}\mathrm{t}\mathrm{o}\mathrm{m}}$,井口加速度${\stackrel{··}{\delta }}_{\mathrm{t}\mathrm{o}\mathrm{p}}$,井底加速度${\stackrel{··}{\delta }}_{\mathrm{b}\mathrm{o}\mathrm{t}\mathrm{t}\mathrm{o}\mathrm{m}}$,其求解流程如图8所示。

将井口与井底加速度、速度、位移分别作差得到井口与井底的位移、速度、加速度滞后量,即

(17)$\left\{\begin{array}{l}\delta_{\mathrm{zhihou}}=\delta_{\mathrm{top}}-\delta_{\mathrm{bottom}} \\\dot{\delta}_{\mathrm{zhihou}}=\dot{\delta}_{\mathrm{top}}-\dot{\delta}_{\mathrm{bottom}} \\\ddot{\delta}_{\mathrm{zhihou}}=\ddot{\delta}_{\mathrm{top}}-\ddot{\delta}_{\mathrm{bottom}}\end{array}\right.$

式(17)中:${\delta }_{\mathrm{z}\mathrm{h}\mathrm{i}\mathrm{h}\mathrm{o}\mathrm{u}}、{\stackrel{·}{\delta }}_{\mathrm{z}\mathrm{h}\mathrm{i}\mathrm{h}\mathrm{o}\mathrm{u}}、{\stackrel{··}{\delta }}_{\mathrm{z}\mathrm{h}\mathrm{i}\mathrm{h}\mathrm{o}\mathrm{u}}$分别为位移滞后量、速度滞后量、加速度滞后量。

5 仿真结果分析

5.1 COMBIN40与BEAM188单元仿真结果对比

由COMBIN40单元来等效抽油杆,式(11)的维度为h×h;由BEAM188单元来等效抽油杆,式(11)的维度为6h×6h。现基于上述计算方法,抽油杆分别采用BEAM188单元和COMBIN40单元等效,施加相同的边界条件,对计算结果进行对比,其数值如表2所示。

由表2可知,抽油杆采用BEAM188单元来等效计算,需要19 h,且在计算过程中结果易发散,而采用COMBIN40单元来等效计算,只需要1 h,计算结果稳定且易收敛,而且两种单元计算结果变化不大,故采用COMBIN40单元等效抽油杆较BEAM188单元,大幅提高了计算效率。

5.2 模型验证

以大庆油田D井为例进行滞后分析,井深1 527.10 m,造斜点250.00 m,最大井斜38.35°,如图9所示;根据图9中选取其中一组参数进行求解,其具体参数如表3所示。为模拟抽油机在举升过程中的往复运动,根据实际施工参数进行边界条件的施加,在井口施加一个已知位移边界,如图10所示,在0~T/2内,模拟抽油杆上提过程,T/2~T内模拟下行过程。从井口到井底的杆柱组配为φ19 mm×59×539.08 m、φ16 mm×88×799.69 m、φ57 mm×1×50.00 m。抽油杆材质为J55钢级,井眼实际轨迹如图11所示,基于上述计算方法,将计算出的载荷从井底到井口进行迭代,从而得到地面示功图,并将计算与实测地面示功图进行了对比,如图12所示,还对实测与计算井口位移进行了对比,如图13所示。

从图12可以看出,仿真计算得到的地面示功图与实测示功图非常接近,各个阶段的相对误差均在±10%以内;其中,上冲程实测和计算最大载荷分别为67.20、64.15 kN,两者误差率为4.54%;下冲程实测和计算载荷分别为32.21、32.13 kN,两者误差率为0.03%。AB段为弹性变形阶段,由于抽油杆在实际施工中开始就存在弹性变形,表现出计算值略大于实测值;BC上行阶段和FG为下行阶段,实测值与计算值出现误差,表明Newmark法在计算时受到摩擦系数改变的影响;CD为卸载阶段,实测值与计算值吻合很好;DE段为排气阶段,泵筒内存在大量气体,抽油杆柱下行过程中压缩空气,导致卸载缓慢;EF为继续卸载阶段。各个阶段的计算结果与实测结果吻合得很好,是因为Newmark法α=0.5,β=0.25时的无条件稳定的特点。从图12还可看出,抽油杆柱上提初始点A和下行终点G不重合,本质是杆柱弹性变形的动态累积与杆管摩阻力的能量损耗不一致的结果,直接体现了井底柱塞的运动滞后于井口。

由图13可以看出,实测与计算井口位移十分接近。以该井的实际施工数据作为边界条件,通过Newmark法和有限元法相结合的方法对该井进行有限元仿真计算,从而把图12中的滞后数据进行量化,得到该井的滞后时间为182.40 s,滞后位移为4.04 m。

5.3 计算参数

依据上述计算方法,对大庆油田的8口井进行仿真计算,8口井对应的最大井斜角依次分别为10.10°、12.99°、30.06°、31.14°、33.98°、36.25°、38.35°、48.24°,分别对不同冲程、冲次、黏度、泵径参数作用下的抽油杆杆柱进行下行滞后仿真求解,其计算参数如表4所示。

5.4 不同生产参数对长冲程抽油杆滞后影响分析

5.4.1 冲程

在8口井计算中,井斜角取38.35°,保持冲次、黏度、泵径不变,冲程分别为10、20、30、40、45、48、49、50 m时,计算5个周期,得到抽油杆的轴向位移、下行速度及加速度,分别如图14~图16所示。

由图14可知,抽油杆在运行过程中,在较小冲程(如 10、20 m)时,井口与井底的位移曲线相位差较小,滞后不明显。随着冲程的增大,滞后位移由冲程10 m时的0.28 m增大到50 m时的23.34 m,即冲程越大滞后位移也越大;随着冲程的增大井口、中间、井底三处位移曲线的“同步性”变差,说明杆柱在运动过程中的弹性变形和摩擦阻力影响加剧,导致能量传递效率下降。

由图15可知,抽油杆柱下行过程中,冲程为10 m时,当井口下行速度达到0.026 1 m/s时,这时井底的速度为0.022 3 m/s,速度相差0.003 8 m/s;冲程为50 m时,井口下行速度达到0.130 9 m/s,这时井底的速度为0.038 7 m/s,速度相差0.092 2 m/s;随着冲程的增大,抽油杆柱下行滞后速度由0.003 8 m/s增大到0.092 2 m/s。亦表明在大斜度长冲程低冲次井中,冲程越大滞后越严重。

由图16可知,也是随着冲程的增大,井口的加速度由0.000 137 m/s²增大到0.000 684 m/s²,井口与井底加速度差值由0.000 018 m/s²增大到0.000 639 m/s²,表明冲程越大滞后越严重。

在8口井计算中,井斜角取38.35°,保持冲次、黏度、泵径不变,冲程分别为10、20、30、40、45、48、49、50 m时,计算5个周期,得到抽油杆井口速度、井底速度及累积摩阻力,如图17所示,其数值如表5所示。

由图17和表5可知,随着冲程的增大,井口最大下行速度随着冲程的增大而增大,其数值由0.026 1 m/s增加到0.130 9 m/s;井底的下行速度数值由0.022 3 m/s增加到0.038 7 m/s,即随着冲程的增大井底和井口的速度都在增大,且累计摩阻力数值也在增大,其数值由0.878 0 kN增加到1.023 0 kN,累计摩阻力越大则抽油杆柱下行滞后越严重。在大斜度长冲程低冲次井的施工中合理降低冲程更能保证采油效率,此结论可为现场施工提供理论指导。

5.4.2 冲次

井斜角取38.35°,保持冲程、黏度、泵径不变,冲次分别为0.020、0.033、0.035、0.050、0.057、0.062、0.065、0.067、0.087 次/min,计算5个周期T得到抽油杆井口速度、井底速度及累积摩阻力,如图18所示,其数值如表6所示。

由图18和表6可知,随着冲次的增加,井口和井底的速度也随之增加,井口速度由0.052 4 m/s增加到0.227 7 m/s,井底速度由0.026 0 m/s增加到0.041 4 m/s,累计摩阻力由0.657 0 kN增加到1.465 0 kN,即随着冲次的不断增加,累计摩阻力也在不断增大,井口与井底的速度差值也在增大,其数值由0.026 4 m/s增加到0.186 3 m/s。综上可知,冲次越大滞后越严重,应合理降低冲次,以减小冲次对抽油杆柱滞后的影响,从而增加采油效率,此结果可为现场施工提供理论指导。

5.4.3 黏度

井斜角取38.35°,保持冲程、冲次、泵径不变,黏度分别为3.8、4.8、5.8、6.8、7.8、8.8、9.8、10.8 mPa·s,计算5个周期T直至完全稳定,得到抽油杆井口速度、井底速度及累积摩阻力,如图19所示,数值如表7所示。

由图19和表7可知,随着黏度的增加,井口最大下行速度基本不变,而井底的速度由0.058 4 m/s减小到0.013 4 m/s,且井口的速度全都要比井底的速度大,原因是抽油杆下行滞后。从图19还可知,随着黏度的增加累计摩阻力也在增大,其数值由0.650 0 kN增加到2.124 1 kN,即随着黏度的不断增加,累计摩阻力增大是导致抽油杆下行滞后的主要原因,黏度越大滞后越严重,在现场施工中应采取相应的措施来降低原油的黏度,如电加热装置等,进而降低原油黏度带来的抽油杆滞后影响。

5.4.4 泵径

井斜角取38.35°,保持冲程、冲次、黏度不变,泵径分别为38、44、57 mm,计算5个周期T直至完全稳定,得到抽油杆的下行速度及下行累计摩阻力结果,其结果如图20所示,数值如表8所示。

由图20和表8可知,随着泵径的增大,井口下行的最大速度基本不变,井底的速度由0.057 5 m/s减小到0.036 7 m/s,原因是随着泵径的逐渐增大,下行的累计摩阻力也在增大,其数值由0.763 kN增加到1.023 kN,即泵径越大,下行滞后越严重。

5.4.5 井斜角

保持冲程、黏度、泵径,冲次不变,井斜角分别为10.10°、12.99°、30.06°、31.14°、33.98°、36.25°、38.35°、48.24°计算5个T。得到抽油杆的下行速度及下行累计摩阻力结果,其结果如图21所示,数值如表9所示。

由图21和表9可知,不同井对应不同的井斜角,随着井斜角的不断增大,井口最大下行速度不变,井底下行速度由0.120 4 m/s减小到0.033 5 m/s,原因是随着井斜角的逐渐增大,下行累计摩阻力也在增大,其数值由0.291 1 kN增加到1.540 2 kN,即井斜角越大,下行滞后越严重,实际施工中应在抽油杆上发生偏磨的位置加扶正器来减少杆管间的摩阻力,从而降低滞后的影响,进而增加采油效率。

6 结论

基于大庆油田某厂提供的8口长冲程井,利用COMBIN40单元模拟抽油杆的运动状态,分析了不同冲程、冲次、采出液黏度、泵径及井斜角参数作用下抽油杆的运动状态。得到以下结论。

(1)采用COMBIN40单元来等效抽油杆,用COMBIN40单元的内置库伦摩擦模型来计算杆管间的摩阻力,并将计算功图与实测功图进行对比,发现两者结果十分接近,即用COMBIN40单元能较好地模拟抽油杆在井下的运动状态,而且用COMBIN40单元来模拟抽油杆,相较于BEAM188单元大幅提高了计算效率。

(2)随着冲程的逐渐增大,滞后位移由冲程10 m时的0.28 m增大到50 m时的23.34 m,滞后速度由10 m时的0.003 8 m/s增大到50 m时的0.092 2 m/s,滞后加速度由10 m时的0.000 018 m/s²增大到50 m时的0.000 64 m/s²,原因是随着冲程的增大杆管累计摩阻力也在增大,其数值由0.878 0 kN增加到1.023 0 kN,抽油杆下行滞后越严重,建议在大斜度长冲程低冲次井中应合理降低冲程,从而增加采油效率,此结果可为现场施工提供理论指导。

(3)对不同冲次下抽油杆运动状态进行仿真计算,得到随着冲次的增加井口和井底的速度也随之增加,冲次越大滞后越严重,而现场提供的8口井均为高含水的低产井,并且泵筒内含有气体,应当合理降低冲次来减小滞后的影响,从而增加采油效率降低能耗。

(4)随着井液黏度的增加,井口最大下行速度基本不变,而井底最大下行速度由0.058 4 m/s减小到0.013 4 m/s,随着黏度的增加累计摩阻力也在增大,其数值由0.650 0 kN增加到2.124 1 kN,随着黏度的不断增加,累计摩阻力增大是导致抽油杆下行滞后的重要原因,在现场施工中应采取相应的措施来降低原油的黏度,如电加热装置来降低原油黏度,进而降低抽油杆滞后对采油的影响。

(5)随着泵径的增大,井口最大下行速度基本不变,井底最大下行速度由0.057 5 m/s减小到0.036 7 m/s,原因是随着泵径的逐渐增大,下行的累计摩阻力也在增大,其数值由0.763 0 kN增加到1.023 0 kN,即泵径越大,下行滞后越严重。

(6)随着井斜角的增大,抽油杆柱下行滞后越严重,井斜角越大表明抽油杆与油管间的接触位置在逐渐增多,实际施工中应在抽油杆上偏磨位置布置扶正器来减少杆管间的摩阻力,从而降低滞后的影响。

参考文献

原文顺序 | 出版日期 | 本文引用

[1]	王宏博, 董世民, 甘庆明, 等. 低产抽油机井动态参数仿真模型与提高系统效率的途径[J]. 石油学报, 2018, 39(11): 1299-1307.

[2]	Wang Hongbo, Dong Shimin, Gan Qingming, et al. Simulation model of dynamic parameters for low-yield pumping units and approaches to enhance system efficiency[J]. Acta Petrolei Sinica, 2018, 39(11): 1299-1307.

[3]	韩修廷, 王德喜, 王研, 等. 利用低摩阻柱塞抽油泵提高泵效及防偏磨技术的应用[J]. 石油学报, 2007, 28(4): 138-141.

[4]	Han Xiuting, Wang Dexi, Wang Yan, et al. Application of low-friction plunger oil pumps to enhance pumping efficiency and prevent eccentric wear[J]. Acta Petrolei Sinica, 2007, 28(4): 138-141.

[5]	Chang X P, Liao J B, Yin H L. Dynamic response analysis of carbon fiber sucker rod column under axial motion[J]. Engineering Research Express, 2025, 7(4): 045565.

[6]	赵岩龙, 方正魁, 邱子瑶, 等. 基于长短时记忆网络的腐蚀工况下抽油杆剩余使用寿命预测[J]. 科学技术与工程, 2021, 21(36): 15429-15433.

[7]	Zhao Yanlong, Fang Zhengkui, Qiu Ziyao, et al. Residual service life prediction of sucker rods under corrosion conditions based on long-short-term memory networks[J]. Science Technology and Engineering, 2021, 21(36): 15429-15433.

[8]	康红兵, 彭振华, 丁雯, 等. 偏磨抽油杆接箍力学性能及承载能力分析[J]. 科学技术与工程, 2023, 23(4): 1526-1533.

[9]	Kang Hongbing, Peng Zhenhua, Ding Wen, et al. Mechanical properties and load-bearing capacity analysis of tapered sucker rod couplings[J]. Science Technology and Engineering, 2023, 23(4): 1526-1533.

[10]	姚飞, 孙玉学, 蔡萌, 等. 大泵油管防喷工艺管柱作业安全性[J]. 科学技术与工程, 2024, 24(3): 1036-1044.

[11]	Yao Fei, Sun Yuxue, Cai Meng, et al. Operational safety of blowout prevention workover strings in large pump oil pipelines[J]. Science Technology and Engineering, 2024, 24(3): 1036-1044.

[12]	李杰, 李屹洋, 苏日古, 等. 新疆油田火驱区块采油井D级杆腐蚀失效分析[J]. 科学技术与工程, 2025, 25(3): 992-998.

[13]	Li Jie, Li Yiyang, Su Rigu, et al. Corrosion failure analysis of grade D rod in production wells in fire drive block of Xinjiang Oilfield[J]. Science Technology and Engineering, 2025, 25(3): 992-998.

[14]	周宇, 冯定, 丁亮亮, 等. 稠油井抽油杆柱黏滞摩擦载荷计算分析[J]. 科学技术与工程, 2019, 19(21): 149-154.

[15]	Zhou Yu, Feng Ding, Ding Liangliang, et al. Calculation and analysis of viscous friction load on sucker rod strings in heavy oil wells[J]. Science Technology and Engineering, 2019, 19(21): 149-154.

[16]	王宏博, 孙秀荣, 李凤鑫, 等. 定向井抽油杆柱纵横扭耦合振动仿真研究[J]. 振动工程学报, 2024, 37(11): 1936-1949.

[17]	Wang Hongbo, Sun Xiurong, Li Fengxin, et al. Simulation study on coupled longitudinal-lateral-torsional vibration of sucker rod strings in directional wells[J]. Journal of Vibration Engineering, 2024, 37(11): 1936-1949.

[18]	Wang X B, Lü L G, Li S, et al. Longitudinal vibration analysis of sucker rod based on a simplified thermo-solid model[J]. Journal of Petroleum Science and Engineering, 2021, 196: 107951.

[19]	Yin J J, Ma H Z. A novel method for predicting the behavior of a sucker rod pumping unit based on the polished rod velocity[J]. Mathematics, 2024, 12(9): 1318.

[20]	付强, 吕瑞典, 陈奇, 等. 超稠油井抽油杆滞后分析与软件设计[J]. 石油机械, 2014, 42(9): 85-88.

[21]	Fu Qiang, Lü Ruidian, Chen Qi, et al. Analysis of sucker rod hysteresis in extra-heavy oil wells and software design[J]. Petroleum Machinery, 2014, 42(9): 85-88.

[22]	孔令路. 抽油机井单周期变冲次技术应用实践与认识[J]. 石油地质与工程, 2025, 39(1): 119-122.

[23]	Kong Linglu. Application practice and insights on single-cycle variable stroke rate technology for pumping unit wells[J]. Petroleum Geology and Engineering, 2025, 39(1): 119-122.

[24]	Yang S Y, Yan Y, Dang N, et al. Dynamic characteristics of high-displacement rod pumping systems in SAGD oil wells[J]. Journal of Petroleum Exploration and Production Technology, 2025, 15(11): 1-17.

[25]	张砚雯, 王旱祥, 于长录, 等. 弯曲载荷下碳/玻复合材料抽油杆失效机制[J]. 中国石油大学学报(自然科学版), 2022, 46(5): 133-140.

[26]	Zhang Yanwen, Wang Hanxiang, Yu Changlu, et al. Failure mechanism of carbon/glass composite sucker rods under bending loads[J]. Journal of China University of Petroleum (Natural Science Edition), 2022, 46(5): 133-140.

[27]	Marino L, Cicirello A. Coulomb friction effect on the forced vibration of damped mass-spring systems[J]. Journal of Sound and Vibration, 2022, 535: 117085.

[28]	王宏博, 董世民, 张洋, 等. 基于气液分离的泵内压力与抽油杆柱纵向振动耦合仿真[J]. 石油学报, 2023, 44(2): 394-404.

[29]	Wang Hongbo, Dong Shimin, Zhang Yang, et al. Simulation of coupling between pump internal pressure and longitudinal vibration of sucker rod string based on gas-liquid separation[J]. Acta Petrolei Sinica, 2023, 44(2): 394-404.

[30]	张强. 小型柔式抽油系统力学分析及应用[D]. 大庆: 大庆石油学院, 2008.

[31]	Zhang Qiang. Mechanical analysis and application of small flexible pumping systems[D]. Daqing: Daqing Petroleum Institute, 2008.

[32]	王新敏, 李义强, 许宏伟. ANSYS结构分析单元与应用[M]. 北京: 人民交通出版社, 2011.

[33]	Wang Xinmin, Li Yiqiang, Xu Hongwei. ANSYS structural analysis elements and applications[M]. Beijing: People’s Communications Press, 2011.

[34]	吴永礼. 计算固体力学方法[M]. 北京: 科学出版社, 2003.

[35]	Wu Yongli. Computational solid mechanics method[M]. Beijing: Science Press, 2003.

[36]	Newmark N M. A method of computation for structural dynamics[J]. Journal of Engineering Mechanics, 1959, 85: 67-94.

基金资助

国家自然科学基金(52374034)

国家自然科学基金(51904075)

黑龙江省自然科学基金(LH2022E018)

PDF (20368KB)

访问

被引

详细

导航

Received	Published	Revised
2025-09-26	2026-08-18	2025-12-18
Online Date
2026-09-04

摘要