考虑分布轴向力的细长杆横向振动与失稳分析1)

图1 抽油杆柱横向弯曲力学模型

2 分布轴向力下杆柱失稳模型及临界载荷计算方法

根据图1力学模型,建立细长杆柱的数学模型为

(1)$\begin{eqnarray} \label{eq1} EI\frac{\partial^{4}y}{\partial x^{4}}+[-p_{0} -q(L-x)]\frac{\partial ^{2}y}{\partial x^{2}}+q\frac{\partial y}{\partial x}=0 \end{eqnarray} $

方程(1)进一步化为

(2)$\begin{eqnarray} \label{eq2} \frac{\partial^{4}Y}{\partial X^{4}}+[-p-(1-X)\beta]\frac{\partial ^{2}Y}{\partial X^{2}}+\beta \frac{\partial Y}{\partial X}=0 \end{eqnarray} $

其中

$\begin{eqnarray*} X=\frac{x}{L},\ \ Y=\frac{y}{L},\ \ p=\frac{p_{0} L^{2}}{EI},\ \ \beta =\frac{qL^{3}}{EI} \end{eqnarray*}$

式中,$L$为抽油杆柱全长,m;$q$为抽油杆柱轴向分布力,N/m;$E$为抽油杆柱材料弹性模量,GPa;$\varepsilon$为抽油杆柱线性分布轴向力比例系数;$I$为抽油杆柱的抗弯惯性矩,m$^{4}$;$p_{0}$为抽油杆柱端部受力(向下为正,受拉;向上为负,受压),N。

边界条件为

(3)$\begin{eqnarray} \label{eq3} \left.\begin{array}{l} \left. {Y\left( {X,t} \right)} \right|_{X=0} =0,\ \ \left. {\dfrac{\partial Y\left( {X,t} \right)}{\partial X}} \right|_{X=0} =0 \\[5mm] \left. {Y\left( {X,t} \right)} \right|_{X=1} =0,\ \ EI\left. {\dfrac{\partial ^{2}Y\left( {x,t} \right)}{\partial X^{2}}} \right|_{X=1} =0 \\ \end{array} \right\} \end{eqnarray} $

式(3)即为含有分布轴向力的杆柱屈曲方程,难于求出精确的解析解。

当$q=0$时,方程(1)化为

(4)$\begin{eqnarray} \label{eq4} EI\frac{\partial^{4}y}{\partial x^{4}}-p_0\frac{\partial^{2}y}{\partial x^2}=0 \end{eqnarray} $

$p_{0}$为负值时,杆柱受压,式(4)与材料力学中的压杆稳定性方程一致。

2.1 有限差分法求解

将杆柱划分为$n$段,得到$n+1$个节点。差分格式为

(5)$\frac{\partial Y}{\partial X}=\frac{Y_{i+1} -Y_{i-1} }{2h} $

(6)$\frac{\partial^{2}Y}{\partial X^{2}}=\frac{Y_{i+1} -2Y_{i} +Y_{i-1} }{h^{2}} $

(7)$\frac{\partial^{4}Y}{\partial X^{4}}=\frac{6Y_{i} -4(Y_{i-1} +Y_{i+1} )+Y_{i-2} +Y_{i+2} }{h^{4}} $

将式(5)$\sim$式(7)代入方程(2)离散化得

(8)$\begin{eqnarray} \label{eq8} &&-p(-h^{2}Y_{i+1} +2h^{2}Y_{i} -h^{2}Y_{i-1} )=\\[2mm]&&\qquad Y_{i+2} +a_{i} Y_{i+1} +b_{i} Y_{i} +c_{i} Y_{i-1} +Y_{i-2} \end{eqnarray} $

式中,$a_{i} ={\beta h^{3}}/{2}-h^{2}(1-X)\beta -4$, $b_{i}=2h^{2}(1-X)\beta +6$, $c_{i} =-{\beta h^{3}}/{2}-h^{2}(1-X)\beta -4$。

由边界条件可知

(9)$\begin{eqnarray} \label{eq9} \left. {\begin{array}{l} Y_{0} =0,\ \ Y_{-1}=Y_{1} \\[1mm] Y_{n} =0,\ \ Y_{n+1} =-Y_{n-1}\\ \end{array}} \right\} \end{eqnarray} $

将式(8)和式(9)展开成矩阵形式得

(10)$\begin{eqnarray} \label{eq10} &&p\left[ {{\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} {-2h^{2}} & {h^{2}} & 0 & \cdots & \cdots & \cdots & \cdots & 0 \\ {h^{2}} & {-2h^{2}} & {h^{2}} & 0 & \cdots & \cdots & \cdots & 0 \\ 0 & {h^{2}} & {-2h^{2}} & {h^{2}} & 0 & \cdots & \cdots & 0 \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0 & \cdots & \cdots & 0 & {h^{2}} & {-2h^{2}} & {h^{2}} & 0 \\ 0 & \cdots & \cdots & \cdots & 0 & {h^{2}} & {-2h^{2}} & {h^{2}} \\ 0 & \cdots & \cdots & \cdots & \cdots & 0 & {h^{2}} & {-2h^{2}} \\ \end{array} }} \right]\left[ {{\begin{array}{*{20}c} {Y_{1} } \\ {Y_{2} } \\ {Y_{3} } \\ \cdots \\ \cdots \\ {Y_{n-3} } \\ {Y_{n-2} } \\ {Y_{n-1} } \\ \end{array} }} \right]=\\&&\qquad \left[ {{\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} {b_{1} +1} & {a_{1} } & 1 & 0 & \cdots & \cdots & \cdots & 0 \\ {c_{2} } & {b_{2} } & {a_{2} } & 1 & 0 & \cdots & \cdots & 0 \\ 1 & {c_{3} } & {b_{3} } & {a_{3} } & 1 & 0 & \cdots & 0 \\ 0 & 1 & {c_{4} } & {b_{4} } & {a_{4} } & 1 & 0 & 0 \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0 & \cdots & 0 & 1 & {c_{n-3} } & {b_{n-3} } & {a_{n-3} } & 1 \\ 0 & \cdots & \cdots & 0 & 1 & {c_{n-2} } & {b_{n-2} } & {a_{n-2} } \\ 0 & \cdots & \cdots & \cdots & 0 & 1 & {c_{n-1} } & {b_{n-1} +1} \\ \end{array} }} \right]\left[ {{\begin{array}{*{20}c} {Y_{1} } \\ {Y_{2} } \\ {Y_{3} } \\ {Y_{4} } \\ \cdots \\ {Y_{n-3} } \\ {Y_{n-2} } \\ {Y_{n-1} } \\ \end{array} }} \right] \end{eqnarray} $

即为

(11)$\begin{eqnarray} p{AY}={BY} \end{eqnarray} $

式(11)两边同乘$A^{-1}$,可化简为$pY=A^{-1}BY$,求出$A^{-1}B$的最小特征值,即可求得$p$,则临界载荷$F_{\rm cr}=-p_{0} =-pEI/L^{2}$。

2.2 伽辽金法

设满足边界条件(3)的形函数为

(12)$\begin{eqnarray} \label{eq12} Y(X)=\sum\limits_{n=1}^n {a_{n} } \varphi_{n} (X) \end{eqnarray} $

其中

(13)$\begin{eqnarray} \left.\begin{array}{l} \varphi_{n} (X)={\rm ch}\lambda_{n} X-\cos \lambda_{n} X-\\ b_{n} ({\rm sh}\lambda_{n} X-\sin \lambda_{n} X) \\ b_{n} =\dfrac{{\rm ch}\lambda_{n} +\cos \lambda_{n} }{{\rm sh}\lambda_{n} +\sin \lambda_{n} },\ \ n=1,2,3\cdots \\ \end{array} \right\} \end{eqnarray} $

式中,

$\lambda_{n} =\beta_{n} L$由超越函数$\cos \lambda_{n} {\rm ch}\lambda_{n} =1$确定,$\lambda_{n} \approx \left( {n+{1}/{4}} \right)\pi $。

同理,得到迦辽金方程为

(14)$\begin{eqnarray} \label{eq13} \int_0^L {L[Y(X)]} \left\{ {\begin{array}{c} \varphi_{1} (X) \\ \varphi_{2} (X) \\ \vdots \\ \end{array}} \right\}{\rm d}X=0 \end{eqnarray} $

整理方程(14),得到关于一待定系数列阵[$a_{n} $]的齐次线性方程组

(15)$\begin{eqnarray} \label{eq15} [L][a_{n} ]=0 \end{eqnarray} $

其中,$L_{nn} =(\lambda_{n} )^{4}I_{1nn} -pI_{2nn} -\beta I_{3nn} $,$L_{mn} =(\lambda_{n} )^{4}I_{1mn} -pI_{2mn} -\beta I_{3mn}+\beta I_{4mn}$。$I_{1nn}$, $I_{2nn}$, $I_{3nn}$,\linebreak $I_{1mn}$, $I_{2mn}$, $I_{3mn}$, $I_{4mn}$为积分常数。

$a_{n} $有非零解的条件为

(16)$\begin{eqnarray} \det L(p)=0\ \ \mbox{或}\ \ \left| {L(p)} \right|=0 \end{eqnarray} $

可求得失稳时的临界载荷,将其代入式(16)可得到振型函数。

3 轴向分布力下的杆柱横向振动数学模型及固有频率计算方法

根据图1力学模型,建立杆柱的横向振动数学模型为

(17)$\begin{eqnarray} \label{eq17} &&\rho A\frac{\partial^{2}y}{\partial t^{2}}+EI\frac{\partial^{4}y}{\partial x^{4}}+q\frac{\partial y}{\partial x}-\\&&\qquad (p_{0} +qL-qx)\frac{\partial ^{2}y}{\partial x^{2}} =0 \end{eqnarray} $

边界条件为

(18)$\begin{eqnarray} \label{eq18} \left.\begin{array}{l} \left. {y\left( {x,t} \right)} \right|_{x=0} =0,\ \ \left. {\dfrac{\partial y\left( {x,t} \right)}{\partial x}} \right|_{x=0} =0 \\[4mm] \left. {y\left( {x,t} \right)} \right|_{x=L} =0,\ \ EI\left. {\dfrac{\partial^{2}y\left( {x,t} \right)}{\partial x^{2}}} \right|_{x=L} =0 \\ \end{array} \right\} \end{eqnarray} $

3.1 积分法

先将$y(x,t)$分离变量,即令$y(x,t)=Z(x)F(t)$,则方程(17)可化为

(19)$\begin{eqnarray} \label{eq19} &&EI\frac{{\rm d}^{4}Z(x)}{{\rm d}x^{4}}+q\frac{{\rm d}Z(x)}{{\rm d}x}-\left( {F_{p} +qL-qx} \right)\frac{{\rm d}^{2}Z(x)}{{\rm d}x^{2}}-\\&&\qquad \omega ^{2}\rho AZ(x)=0 \end{eqnarray} $

为便于计算,式(19)中参数进行变换,则

$\begin{eqnarray*} \left.\begin{array}{l} X=\dfrac{x}{L},\ \ Y(X)=\dfrac{Z\left( x \right)}{L},\ \ \varOmega^{2}=\dfrac{\rho A\omega^{2}L^{4}}{EI}\\[3mm] \alpha =\dfrac{qL^{3}}{EI},\ \ r^{2}=\dfrac{p_{0} L^{2}}{EI} \end{array}\right. \end{eqnarray*}$

对式(19)进行四次积分,并整理得

(20)$\begin{eqnarray} \label{eq20} &&Y+\int_0^X {k_{1} (X,s)z(s){\rm d}s} +\\&&\qquad \int_0^1 {k_{2} (X,s)Y(s){\rm d}s} =0 \end{eqnarray} $

式中

(21)$\begin{eqnarray} \left.\begin{array}{l} k_{1} (X,s)=-\left( {X-s} \right)\bigg[\dfrac{\alpha }{2}\left( {X-3s+2} \right)+\\[3mm]\qquad \dfrac{\varOmega^{2}}{6}\left( {X-s} \right)^{2}+r^{2}\bigg] \\[3mm] k_{2} (X,s)=\dfrac{X^{3}}{4}\bigg\{\alpha +\varOmega^{2}\left( {1-s} \right)-\\[3mm]\qquad 2\bigg[\dfrac{\alpha }{2}\left( {3-3s} \right)+\dfrac{\varOmega^{2}}{6}\left( {1-s} \right)^{2}+r^{2}\bigg]\cdot\\[3mm]\qquad\left( {1-s} \right)\bigg\}+ \dfrac{X^{2}}{4}\bigg\{6\bigg[\dfrac{\alpha }{2}\left( {3-3s} \right)+\\[3mm]\qquad\dfrac{\varOmega ^{2}}{6}\left( {1-s} \right)^{2}+r^{2}\bigg]\left( {1-s} \right)- \alpha -\\[3mm]\qquad\varOmega ^{2}\left( {1-s} \right)\bigg\} \\[3mm] \end{array} \right\} \end{eqnarray} $

设$z(\xi )$为多个多项式之和,采用级数解来逼近其真实解,则

(22)$\begin{eqnarray} \label{eq22} Y\left( X \right)=\sum\limits_{n=0}^{N+2} {c_{n} X^{n+2}} ,\ \ 0\leqslant X\leqslant 1 \end{eqnarray} $

将式(22)代入式(20)得

(23)$\begin{eqnarray} \label{eq23} &&\sum\limits_{n=0}^{N+2} {\left( {b_{mn} +k_{1mn} +k_{2mn} }\right)c_{n} } =0\\&&\qquad (m=0,1,2, \cdots, N) \end{eqnarray} $

式中

$\begin{eqnarray*} \left.\begin{array}{l} b_{mn} =\dfrac{1}{m+n+3}\\ k_{1mn} =\int_0^1 X^{m}\int_0^X k_{1} \left( {X,s} \right)s^{n+2}{\rm d}s{\rm d}X \\ k_{2mn} =\int_0^1 {X^{m}\int_0^1 {k_{2} \left( {X,s} \right)s^{n+2}{\rm d}s{\rm d}X} } \end{array}\right. \end{eqnarray*}$

若$c_{n}$有非零解,则

(24)$\begin{eqnarray} \label{eq24} &&\det \left( {b_{mn} +k_{1mn} +k_{2mn} } \right)=0\\&&\qquad (m,n=0,1,2, \cdots, N+2) \end{eqnarray} $

由式(24)求解出$\varOmega $,即可进一步得到杆柱横向振动固有频率$\omega $。特殊地,当$q=0$时,$p_{0}=0$,方程(17)化为

(25)$\begin{eqnarray} \label{eq25} \rho A\frac{\partial^{2}y}{\partial t^{2}}+EI\frac{\partial^{4}y}{\partial x^{4}}=0 \end{eqnarray} $

式(25)与教材^[2]中梁的横向自由振动方程是相同的,此方程有解析解,亦可采用此方法求解。

3.2 伽辽金法

令$\varOmega t=\varOmega \sqrt {{\rho AL^{4}}/({EI})} T=\omega T$,$\omega $为抽油杆柱的振动频率。其他参数同方程(3)。对方程(17)进行参数变换得

(26)$\begin{eqnarray} \label{eq26} &&\frac{\partial^{4}Y}{\partial X^{4}}+[(-p-\beta(1-X)]\frac{\partial^{2}Y}{\partial X^{2}}+\\&&\qquad \beta \frac{\partial y}{\partial x}+\frac{\partial^{2}Y}{\partial T^{2}}=0 \end{eqnarray} $

令

(27)$\begin{eqnarray} \label{eq27} Y(X,T)=Y(X)\cos (\omega T) \end{eqnarray} $

将式(27)代入方程(17),得振型方程为

(28)$\begin{eqnarray} \label{eq28} &&\frac{\partial^{4}Y}{\partial X^{4}}+[(-p-\beta(1-X)]\frac{\partial^{2}Y}{\partial X^{2}}+\\&&\qquad \beta \frac{\partial Y}{\partial X}-\omega^{2}Y=0 \end{eqnarray} $

设满足边界条件的形函数为梁自由振动函数

(29)$\begin{eqnarray} \label{eq29} Y(X)=\sum\limits_{n=1}^n {a_{n} } \varphi_{n} (X) \end{eqnarray} $

$\varphi_{n} (X)$同式(13)。将式(28)和式(29)代入式(27),得到迦辽金方程为

(30)$\begin{eqnarray} \label{eq30} \int_0^L {L[Y(X)]} \left\{ {\begin{array}{c} \varphi_{1} (X) \\ \varphi_{2} (X) \\ \vdots \\ \end{array}} \right\}{\rm d}X=0 \end{eqnarray} $

整理方程(30),得到关于一待定系数列阵[$a_{n} $]的齐次线性方程组为

(31)$\begin{eqnarray} \label{eq31} [L][a_{n} ]=0 \end{eqnarray} $

其中

$\begin{eqnarray*} \left.\begin{array}{l} L_{nn} =\big[(\lambda_{n})^{4}-\omega ^{2} \big]I_{1nn}- pI_{2nn}-\beta I_{3nn}\\ L_{mn} =(\lambda_{n} )^{4}I_{1mn}-pI_{2mn}-\beta I_{3mn}+\beta I_{4mn} \end{array}\right. \end{eqnarray*}$

$I_{1nn}$, $I_{2nn}$, $I_{3nn}$, $I_{1mn}$, $I_{2mn}$,$I_{3mn}$, $I_{4mn}$为积分常数。$a_{n} $有非零解的条件是

(32)$\begin{eqnarray} \det L(\omega )=0\ \ \mbox{或}\ \ \left| {L(\omega )} \right|=0 \end{eqnarray} $

可求得各阶固有频率,将其代入式(32)可得到振型函数。

4 两种模型数值分析

4.1 杆柱失稳临界载荷

4.1.1 准确性测试

当分布轴向力$q=0$时,数学方程(1)简化为《材料力学》教材^[1]中的模型,而数学方程(17)简化为《机械系统动力学》教材^[2]中的模型。其他参数为$L=5$ m,$d=20$ mm,$E=209$ GPa时,材料力学欧拉杆临界载荷公式为

(33)$\begin{eqnarray} \label{eq33} F_{\rm cr} =\frac{\pi^{2}EI}{(0.7L)^2} \end{eqnarray} $

采用式(33)欧拉法、2.1节差分法、2.2节伽辽金法分别求解相同参数下杆柱的临界载荷为1322.5 N, 1325.7 N, 1321.2 N,可知,差分法和伽辽金法的结果与公式(33)求解结果误差微小(误差在0.25%以内),说明差分法和伽辽金法求解欧拉杆临界载荷精度较高。

4.1.2 考虑分布轴向力失稳载荷

取基本参数:$E=209$ GPa,$d=20$ mm。差分法和伽辽金法得出不同杆长下的杆柱临界载荷规律如图2和图3所示。由两图可知,在杆长为50m以内时,差分法和伽辽金法得到的临界载荷规律相同,即随杆长的增加而减小;但杆长为50 m以上时,伽辽金法所得杆柱临界载荷会随杆长的增加而增加,这与差分法所得规律截然相反。由此可知,两种方法在计算考虑分布轴向力较长杆柱的临界载荷时有一种方法得到的结果是不准确且不适用的。而伽辽金方法实质是一种近似计算方法,在杆长超越极限时结果会失真,因而不再适用。这也说明近似方法会有它本身的适用范围,而伽辽金法不适用于考虑分布轴向力且杆长超越一定长度的杆柱的临界载荷的计算。

图2

图2 差分法所得临界载荷变化

图3

图3 伽辽金法所得临界载荷变化

4.2 杆柱横向振动固有频率和振型

4.2.1 准确性测试

情况一:当基本参数为$L=5$ m,$d=20$ mm{,} $E=209$ GPa,$q=0$ N/m,且底部载荷$p_{0}=0$时,杆柱则做横向自由振动,其固有频率可表示^[2]为

(34)$\begin{eqnarray} \label{eq34} \omega_{r} =\frac{(r+1/4)^{2}\pi^{2}}{L^{2}}\sqrt {\frac{EI}{\rho A}}\ \ (r=1,2,3 \cdots) \end{eqnarray} $

采用式(34)、3.1节积分法和3.2节伽辽金法求解该相同参数下杆柱的固有频率,其结果对比如表1所示。由表1可知,前六阶固有频率所得结果误差很小(0.15%以内)。因此,不考虑端部载荷和分布轴向力时,积分法和伽辽金法求解杆柱横向振动固有频率和振型精度都很高。

表1 未考虑分布轴向力计算固有频率(单位:HZ)

情况二:在情况一参数基础上,只变动$p_{0}$,令$p_{0}=2000$N,得到杆柱横向振动的固有频率,如表2所示。由表2可得,底端杆柱受拉的情况下,积分法和伽辽金法求解杆柱固有频率时误差甚小(0.02%以内),说明两种方法均适应于杆柱受拉情况杆柱固有频率计算。

表2 考虑分布轴向力计算固有频率(单位:HZ)

4.2.2 考虑分布轴向力和杆长变化的横向振动

情况三:固有频率随轴向分布力和杆长的变化规律。假设基本参数为$d=20$ mm,$E=209$ GPa,$p_{0}=500$N,改变轴向分布力$q$和$L$,得到杆柱固有频率随分布轴向力$q$和$L$的关系曲线,如图4所示。由图4可知,固有频率随分布轴向力的增大而增大,而随杆柱长度的增长而减小。$q$变化范围在0$\sim$20 N/m时,50 m内杆柱长度固有频率最大误差率为26.96%,100$\sim$2000 m内杆柱长度固有频率最大误差率为390.6%。由此说明,在杆柱长度较短时,分布轴向力对固有频率的影响不大,杆柱较长时,分布轴向力对固有频率影响较大。杆柱较长时,分布轴向力不宜忽略。

图4

图4 固有频率随杆长$L$和分布轴向力$q$的变化曲线

情况四:固有频率随底端载荷变化的规律。考虑到分布轴向力在杆柱较长时不宜忽略,取基本参数$q=20$ N/m,$L=1000$ m,$d=25$ mm,$E=209$ GPa,得到不同载荷下杆柱的固有频率变化曲线,如图5所示(图中$p_{0}>0$为杆柱受拉,反之受压)。由图可知,随着底部受拉载荷的减小,固有频率逐渐下降,底端载荷为0后继续减小载荷(即反向加载,受压载荷逐渐增大),受压载荷越过一定值后,固有频率变为虚数。说明在底部载荷逐渐变化的过程中,有一载荷对应着固有频率为0的情况。令一阶固有频率为0,由2.1节求得静力学下临界载荷即$p_{0}=-169.8$N,失稳变形曲线如图6所示。由3.1节求得动力学下横向振动的变形曲线同样如图6所示。由图可知,两变形几乎完全重合,杆柱失稳时最大变形均出现在杆柱底部附近。说明此时静力学下考虑分布轴向力$q=20$N/m杆柱失稳的变形即为动力学下横向振动变形,而$p_{0}$则为失稳时的临界载荷。由此可知,考虑分布轴向力的杆柱在非零变形时的固有频率为0,即一阶固有频率为0时,杆柱达到失稳的必要条件。特殊的,当分布轴向力为0且杆柱屈曲时,其固有频率对应为0,该结论与文献[8]一致。

图5

图5 固有频率随底端载荷的变化曲线

图6

图6 杆柱变形规律

图6同时给出了不考虑分布轴向力(即$q=0$)时杆柱的静力失稳变形曲线,与$q=20$N/m时失稳变形曲线对比可知,两者的最大变形位置(即突变处)明显不同,即杆柱最先达到失稳的位置点不同。不考虑分布轴向力的杆柱最先在靠近中间位置失稳,而考虑分布轴向力的杆柱最先在靠近端部位置处失稳(距离底端约为8.5m处),这也说明了分布轴向力对杆柱最先失稳位置有显著的影响,即超长杆柱计算失稳问题须考虑分布轴向力。

5 工程实例求解

5.1 考虑分布轴向力的螺杆泵杆柱固有频率和振型

数学模型方程(17)和方程(18)可应用在求解螺杆泵杆柱的固有频率上,其参数如表3所示。计算的前五阶固有频率如表4,其对应的振型函数如图7所示。

表3 螺杆泵基本参数表

表4 考虑分布轴向力的固有频率结果(单位:HZ)

图7

图7 杆柱各阶振型函数曲线

由表4可知,当细长杆柱底端受拉时($p_{0}>0$),积分法和伽辽金方法得到的固有频率误差很小(0.5%以内),可以忽略不计。图7给出了积分法和伽辽金方法求解前四阶杆柱的振型函数曲线,经对比可知,两种方法用于求解杆柱固有频率和振型函数,误差很小,精度较高。

5.2 有杆抽油系统杆柱失稳情况数值分析

有杆抽油系统抽油杆柱的基本参数如表5所示。实际油井浅则几百米深则几千米,根据前面4.1节的结果分析,伽辽金方法不适用于求解该长度杆柱的临界载荷,只能采用有限差分法。其不同杆长下临界载荷情况如表6和图8所示,其对应失稳变形情况如图9所示。

表5 基本参数

表6 失稳临界载荷

图8

图8 失稳载荷随杆长变化曲线

图9

图9 杆柱不同杆长下的失稳变形曲线

由表6和图8可知,对于杆柱长度300$\sim$2000 m的杆柱,考虑分布轴向力和不考虑分布轴向力,其失稳临界载荷相差甚大,因此油田抽油杆柱求解偏磨临界载荷时,分布轴向力不宜忽略不计。

由图8进一步可知,分布轴向力$q=20$ N/m,杆柱长度为200$\sim$2000 m时,其临界载荷随杆长的增加而逐渐减小,直至趋于一稳定值,变化率区间约为0%$\sim$13.6%。说明200m的杆柱代替2000 m杆柱计算临界载荷时,最大误差为13.6%。由此给出工程建议,在精度要求不太高的情况下,考虑分布轴向力的有杆抽油系统,求解杆管偏磨临界载荷时可以以200m左右的杆柱代替几千米的情况进行计算。

由图9可知,不同井深下,杆柱的失稳变形曲线的规律是近似的,即靠近底端处产生最大失稳变形,这与不考虑分布轴向力时在杆柱的中间位置处产生最大失稳变形有显著不同。因此,抽油杆柱的杆管偏磨严重区域会出现在较深且靠近底端部位,这与文献[16,17]结论是一致的。

6 结论

本文建立了考虑分布轴向力的细长杆失稳和横向振动的力学、数学模型,通过数值分析结果对比,得出以下结论:

(1)考虑分布轴向力的细长杆柱失稳的必要条件是其横向振动的固有频率为0。该结论是对先前教材的知识点的延伸,同时也适用于不考虑分布轴向力的杆柱失稳情况。

(2)计算考虑分布轴向力的杆柱受压失稳载荷时,伽辽金近似计算法在求解较长杆柱时超出允许精度范围导致计算结果失真;差分法则适用于求解较长杆柱的失稳载荷,且精度较高。对于端部受拉的杆柱,数值法和伽辽金法计算杆柱固有频率和振型函数的方法都是适用的。

(3)对于抽油杆柱等超长杆,计算失稳临界载荷时考虑分布轴向力是必要的,反之会产生较大的误差。当考虑分布轴向力的超长杆柱(如200m以上的抽油杆柱),其失稳临界载荷随杆长的变化逐渐缓慢减小。精度要求不高时,可不必取杆柱全长,取一部分杆柱计算临界载荷即可。

(4)本文杆柱与《材料力学》中欧拉杆柱失稳突变处的位置明显不同,这是因为本文模型在欧拉杆的基础上考虑了分布轴向力因素所致。本文对《材料力学》的杆柱静力失稳和《机械系统动力学》的梁的横向振动知识进一步拓展和延伸,并提供了工程应用实例。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

孙训方, 方孝淑, 关来泰. 材料力学. 北京: 高等教育出版社, 2019

Sun

Xunfang

, Fang

Xiaoshu

, Guan

Laitai

Mechanics of Materials. Beijing: Higher Education Press, 2019 (in Chinese)

[2]

贺利乐

. 机械系统动力学. 北京: 国防工业出版社, 2014

[本文引用: 4]

Lile

. Mechanical System Dynamics. Beijing: National Defense Industry Press, 2014 (in Chinese)

[本文引用: 4]

[3]

王平, 白象忠, 刘立静

等.

横向磁场中细长压杆的分岔特性

力学与实践, 2010,32(2):103-107

URL [本文引用: 1]

在磁弹性非线性运动方程、物理方程、电动力学方程及洛仑兹力表达式的基础上，应用Lagrange描述法建立了横向磁场中两端铰支受压细长杆的非线性磁弹性动力学模型. 通过对该模型的简化，分别讨论了静力学模型、线性动力学模型和含三次非线性项的动力学模型的分岔特性. 最后通过数值计算，给出了横向磁场中受压细长杆的失稳临界载荷与相关参量之间的关系曲线，并对计算结果及其变化规律进行了分析讨论.

Wang

Ping

, Bai

Xiangzhong

, Liu

Liling

, et al.

Bifurcation of an Euler's pole in a transverse magnetic field

Mechanics in Engineering, 2010,32(2) : 103-107 (in Chinese)

[4]

侯祥林, 范炜, 贾连光.

变截面压杆临界载荷的迭代算法

哈尔滨工业大学学报, 2011,43(S1):237-240

Hou

Xianglin

, Fan

Wei

, Jia

Lianguang

, et al.

Iterative algorithm of the critical load about variable cross-section bar

Journal of Harbin Institute of Technology, 2011,43(S1):237-240 (in Chinese)

[5]

董冠文, 李宗义, 赵彦军

等.

压杆稳定临界力欧拉公式统一推导

武汉工程大学学报, 2012,34(12):71-74

Dong

Guanwen

, Li

Zongyi

, Zhao

Yanjun

, et al.

Unified deduction of pressure lever stability critical force Euler formula

Journal of Wuhan Institute of Technology, 2012,34(12):71-74 (in Chinese)

[6]

武秀根, 郑百林, 贺鹏飞.

限制失稳杆的后屈曲分析

同济大学学报(自然科学版), 2009,37(1):30-34, 114

Xiugen

, Zheng

Bailin

, He

Pengfei

Study on post-buckling of constrained bar

Journal of Tongji University (Natural Science), 2009,37(1):30-34, 114 (in Chinese)

[7]

丁虎, 邵志华.

弹性结构横向弯曲振动传递率

力学与实践, 2019,41(6):661-664

Ding

, Shao

Zhihua

Transmissibility of transverse bending vibration of elastic structures

Mechanics in Engineering, 2019,41(6):661-664 (in Chinese)

[8]

腾兆春, 李世荣.

轴向力作用下非对称支承 Euler-Bernoulli 梁在过屈曲前后的横向自由振动

振动与冲击, 2004,23(4):88-91

Teng

Zhaochun

, Li

Shirong

Transverse free vibration in pre-buckling and post-buckling of a unsymmetric supports Euler-Bernoulli Beam under axially compressive force

Journal of Vibration and Shock, 2004,23(4):88-91 (in Chinese)

[9]

赵伟东, 杨亚平, 李秀莲

等.

简谐干扰力作用的弹性基础上压杆的横向振动与屈曲

应用力学学报, 2015,32(4):642-646, 708

Zhao

Weidong

, Yang

Yaping

, Li

Xiulian

, et al.

The transverse vibration and buckling of a compressive bar on elastic foundation subjected to harmonic excitation

Chinese Journal of Applied Mechanics, 2015,32(4):642-646, 708 (in Chinese)

[10]

田耀宗, 蹇开林.

轴向运动梁的横向振动分析

应用数学和力学, 2019,40(10):1081-1088

Tian

Yaozong

, Jian

Kailin

Lateral vibration analysis of axially moving beams

Applied Mathematics and Mechanics, 2019,40(10):1081-1088 (in Chinese)

[11]

鲍四元, 曹津瑞, 周静.

任意弹性边界下非局部梁的横向振动特性研究

振动工程学报, 2020,33(2):276-284

Bao

Siyuan

, Cao

Jinrui

, Zhou

Jing

Transverse vibration characteristics of nonlocal beams with arbitrary elastic boundary conditions

Journal of Vibration Engineering, 2020,33(2):276-284 (in Chinese)

[12]

Shin

, Yun

Transverse vibration of a uniform Euler-Bernoulli beam under varying axial force using differential transformation method

Journal of Mechanical Science & Technology, 2006,20:191-196

[13]

Boiangiu

, Costin

, Untaroiu

A transfer matrix method for free vibration analysis of Euler-Bernoulli beams with variable cross section

Journal of Vibration & Control, 2016,22(11):2591-2602

[14]

孙秀荣, 董世民, 王宏博

等.

整体抽油杆柱在油管内空间屈曲的多段式仿真模型对比

吉林大学学报(工学版), 2018,48(4):1124-1132

Sun

Xiurong

, Dong

Shimin

, Wang

Hongbo

, et al.

Multistage simulation models comparison of the entire sucker rod with spatial buckling in the tubing

Journal of Jilin University (Engineering and Technology Edition), 2018,48(4):1124-1132 (in Chinese)

[15]

张雷, 孙秀荣.

铅直圆筒内细长杆柱在屈曲位移激励下的横向振动仿真

力学与实践, 2019,41(5):526-533

Zhang

Lei

, Sun

Xiurong

Transverse vibration simulation of slender rod string in the vertical cylinder under buckling displacement excitation

Mechanics in Engineering, 2019,41(5):526-533 (in Chinese)

[16]

孙秀荣, 江少波, 黄有泉

等.

直井抽油杆柱在屈曲位移激励下横向振动的仿真研究

振动工程学报, 2018,31(5):854-861

Sun

Xiurong

, Jiang

Shaobo

, Huang

Youquan

, et al.

The simulation research of transverse vibration of sucker rod string under buckling deformation excitation in vertical wells

Journal of Vibration Engineering, 2018,31(5):854-861 (in Chinese)