0 引言
近年来，随着高强度钢材的大量使用和自动单缸插销技术的成熟，多节伸缩臂架技术广泛应用到大型工程起重机上。随着钢结构高度的不断增大，钢材的消耗量也在不断增加，梁柱的截面沿轴向发生变化是一种经济型承担载荷的方式。箱形阶梯柱模型以其长度易于伸缩变化，承载力大和材料利用率高等优点,广泛用于工程起重机、高空救援车等各类机械设备中，其稳定性为学者研究的重点[1-2]。TIMOSHENKO等[3]对阶梯柱模型进行了研究和分析，并使用能量法对高阶阶梯柱的临界力进行了计算。王欣等[4]借助预设近似挠度曲线，使用能量法和李兹法对阶梯柱进行了研究，但该方法相当于引入了附加约束，特别是对于四阶以上的阶梯柱会产生较大的误差。文献[5-7]使用传递矩阵法对阶梯柱的临界力进行了研究，但该方法目前也只能用于阶数较少的阶梯柱，在阶梯柱阶数较大情况下使用该方法也会出现较大的误差。随着有限元理论的发展和完善，陆念力等[8-10]使用精确有限元法对阶梯柱的临界力进行了研究。
我国现行的国家标准GB/T3811—2008[11]将文献[8-10]的方法作为规范的阶梯柱稳定性分析精确计算方法，但这种方法对于多阶阶梯柱来说，其刚度矩阵相当庞大，特征方程复杂，对于常见支撑形式的阶梯柱稳定性精确数值解，常通过试凑法获得，计算量巨大，给实际应用带来了一定的困难。在现行的设计规范中，使用图表来表示2～5阶的阶梯柱部分特殊组合情况下的长度系数，但当阶数超过5阶时便无法应用，并且在非特定组合时使用线性拟合的方法，没有给出线性拟合的误差。德国Liebherr LTM 11200、徐工集团QAY1200、三一重工的SANY SAC12000、中联重科的QAY2000等机型的伸缩臂架阶数都达到了8，此时现有的规范就无法满足实际工程的需要,因此需要一种方法，可以较方便且迅速地求解5阶及以上阶梯柱临界力的方法。本文基于纵横弯曲理论建立了n阶伸缩臂阶梯柱的微分方程组，并利用数学归纳法推导出n阶阶梯柱压杆稳定性的递推公式。
1 阶梯柱结构变形及其微分方程
图1 变截面阶梯柱模型及受力简图
Fig.1 Mechanical model and force diagram of a multi-stepped column
针对n阶阶梯柱模型，基于纵横弯曲理论可建立各节伸缩臂挠曲微分方程。如图1所示，n阶阶梯柱的总长为ln，P为伸缩臂顶端承受的轴向力且P的方向保持不变，δ为伸缩臂顶端的侧向位移。假设轴向力和弯矩全部由伸缩臂承受[12-13]，那么第i阶压杆模型的受力和变形如图2所示。
图2 第i阶压杆模型受力及位移图
Fig.2 Force and displacement diagram of thei-th section
根据模型可列出以下方程：
(1)
i=1,2,…,n
式中，Ii为第i节伸缩臂的截面惯性矩，m4；li为第i节伸缩臂顶部到吊臂根部的长度，m；yi为第i节伸缩臂在平面内偏离x轴的侧向位移，m；y″2为yi的二阶导数；E为弹性模量，Pa。
式(1)可统一表示为
(2)
(3)
式(2)的通解为
yi=Aisinkix+Bicoskix+δ
(4)
由挠度位移边界条件可解出各积分常数之间的关系为
(5)
A1=0 B1=-δ
(6)
(7)
(8)
记
(9)
可得到未知系数An和Bn的表达式
(10)
由伸缩臂顶部边界条件: x=ln时yn=δ，得
Ansinknln+Bncosknln=0
(11)
将式(10)代入式(11)中得到伸缩臂失稳特征方程为
(12)
对于一个特定起重机箱形伸缩臂，将所有已知条件代入失稳特征方程可知，式(12)是以P为未知量的非线性方程，即失稳特征方程为超越方程，解此非线性方程即可求得结构失稳临界力P。
2 超越方程的递推规律
当n=1时，根据式(12)可以得到
cosk1l1=0
(13)
当n=2时，可以得到
(14)
当n=3时，可以得到
(15)
根据数学归纳法，可以证明n阶阶梯柱的失稳特征方程如下：
sinknlnC(n)+cosknlnD(n)=0
(16)
C(n)=(sinkn-1ln-1sinknln-1+
(coskn-1ln-1sinknln-1-
(17)
D(n)=(coskn-1ln-1cosknln-1+
(sinkn-1ln-1cosknln-1-
(18)
D(1)=1 C(1)=0
(19)
3 超越方程的数值解法
3.1 超越方程的非线性特征
已知某二阶阶梯柱, l1=15.4 m，l2=29.9 m，I1=43.26×10-3 m4，I2=26.09×10-3 m4，利用式(14)可以得到二阶阶梯柱所构成函数的曲线如图3所示。此函数是一非线性函数且不具有规则的周期性，满足失稳特征方程f(k1,k2)=0的点较多,所以必须通过增加方程的方法来确定临界点。将式(3)进行变形可得增加的方程如下:
(20)
图3 二阶阶梯柱的特征方程多项式与k1函数关系
Fig.3 Relationship between k1 and buckin g characteristic equation of the stepped column wit h two sections
对于此二阶阶梯柱，数值求解方程组为
数值求解结果如下：k1=0.048 0，k2=0.061 7；得到临界力P=2.049 1×107N，长度系数μ2=1.060 12。
3.2 n阶超越方程组的非线性数值解法
当阶梯柱的阶数为n时，超越方程组可以通过增加其约束限制方程建立方程组来求解。所建立的n阶超越方程组如下：
(21)
到目前为止，这种非线性超越方程还不存在解析解，只有数值解。数值解常使用的方法包括欧拉法、龙格库塔法、Gauss-Newton法[14]等，但这些方法可能由于矩阵的奇异而无解。本文使用Levenberg-Marquardt算法[14]来进行求解,这种算法具有梯度法和牛顿法的优点并较原来的梯度下降法求解速度提高了几十倍甚至上百倍，并且由式(21)所构建的方程组总能得到最优数值解。
通过式(21)可以求解出每节阶梯柱的刚度，进而可以求出整个阶梯柱模型的临界力P，此临界力与相同约束条件下与阶梯柱长度相等的基本臂截面所构成均等截面柱的临界力对比，就可以求出n阶阶梯柱的长度系数μ2。计算公式如下：
(22)
3.3 二阶阶梯柱的非线性及其插值
由图4可以看出，使用本文数值方法计算的二阶阶梯柱长度系数与使用ANSYS 17.0计算出的长度系数差值非常小，最大相对误差为1.7×10-4，由此可知n阶压杆稳定性计算的递推公式精度较高。此外，随着第二阶柱的惯性矩的增大，长度系数μ2减小，即阶梯柱的临界力逐渐增大；随着第一阶臂架长度l1与全长臂架长度l之比α1增大即第二阶柱的相对长度减小，长度系数μ2非线性减小。因此在使用GB/T3811—2008中的方法查表得到μ2值时，如果选用了表中没有的βi，那么需使用数值法进行线性插值才可以得到更为准确的结果；另外α1也会影响最终的长度系数值。
图4 二阶阶梯柱的长度系数μ2值
Fig.4 Length coefficient μ2 of stepped colum n with two sections
3.4 多阶阶梯柱的长度系数
表1～表3所示为不同方法计算的三阶、四阶、五阶阶梯柱长度系数μ2的比较。使用数值法对某8阶阶梯柱进行临界力计算，其参数和结果如表4所示。最终可以得到临界力P=1.3049×106 N ,长度系数μ2=2.680。
表1 三阶阶梯柱计算长度系数及比较
Tab.1 Calculated length factor and comparison of3 stepped column
β2β3数值法GB/T3811—2008ANSYS 17.0本文数值法与ANSYS 17.0方法相对误差1.31.31.052 961.0531.053-3.798 70×10-51.31.61.061 761.061 7601.31.91.070 731.070 7301.32.21.079 911.079 99.260 12×10-61.32.51.089 291.0891.089 271.836 09×10-51.61.31.099 021.0991.099 019.099 10×10-51.61.61.109 911.109 891.801 98×10-51.61.91.121 061.121 041.784 06×10-51.62.21.132 471.132 433.532 23×10-51.62.51.144 111.1441.144 073.496 29×10-51.91.31.144 301.1441.144 272.621 76×10-51.91.61.157 231.157 202.592 46×10-51.91.91.170 491.170 444.271 90×10-51.92.21.184 031.183 984.223 04×10-51.92.51.197 841.1981.197 785.009 27×10-52.21.31.188 681.1891.188 652.523 87×10-52.21.61.203 591.203 544.154 41×10-52.21.91.218 851.218 794.922 92×10-52.22.21.234 431.234 356.481 14×10-52.22.51.250 311.2501.250 217.998 66×10-52.51.31.232 121.2321.232 064.869 89×10-52.51.61.248 901.248 844.804 46×10-052.51.91.266 071.265 996.319 17×10-52.52.21.283 591.283 497.791 26×10-52.52.51.301 421.3011.301 309.221 55×10-5
表2 四阶阶梯柱计算长度系数及比较
Tab.2 Calculated length factor and comparison of4 stepped column
β2β3β4数值法GB/T3811—2008ANSYS 17.0本文数值法与ANSYS 17.0方法相对误差1.61.31.31.147 221.1471.147 31-7.844 44×10-51.91.31.31.206 411.2071.206 54-1.077 46×10-42.21.31.31.263 571.2641.263 75-1.424 33×10-42.51.31.31.318 771.3191.319 02-1.895 35×10-41.31.61.31.112 581.1131.112 66-7.189 98×10-51.31.91.31.139 531.1401.139 62-7.897 37×10-51.32.21.31.166 681.1671.166 79-9.427 57×10-51.32.51.31.193 931.1941.194 07-1.172 46×10-41.31.31.61.090 661.090 73-6.417 72×10-51.31.31.91.095 451.095 53-7.302 40×10-51.31.32.21.100 341.100 43-8.178 62×10-51.31.32.51.105 331.1051.105 43-9.046 25×10-5
表3 五阶阶梯柱计算长度系数及比较
Tab.3 Calculated length factor and comparison of5 stepped column
β2β3β4β5数值法GB/T3811—2008ANSYS 17.0本文数值法与ANSYS 17.0方法相对误差1.31.61.31.31.205 871.2061.205 868.292 84×10-61.31.91.31.31.258 501.2591.258 5001.32.21.31.31.310 011.3101.310 0101.32.51.31.31.360 301.3601.360 3001.61.31.31.31.239 701.2401.239 71-8.066 40×10-61.91.31.31.31.322 121.3221.322 092.269 13×10-62.21.31.31.31.400 101.4001.400 072.142 75×10-52.51.31.31.31.474 201.4741.474 21-6.783 29×10-6
表4 某八阶阶梯柱临界力计算参数及结果
Tab.4 Parameters and results of a 8 steppe d column critical force
阶数li(m)Ii(m4)ki115.50043.260 0×10-30.012 1227.12533.276 9×10-30.013 8338.75025.597 6×10-30.015 7450.37519.690 0×10-30.017 9562.00015.146 5×10-30.020 5673.62511.651 2×10-30.023 3785.2508.962 44×10-30.026 6896.8756.894 19×10-30.030 3
八阶阶梯柱计算长度系数及比较见表5，其中μ为由多级阶梯柱侧向刚度分析实用算式得到的计算长度系数[11],μa为使用精确有限元法得到的等效构件的计算长度系数[9]，Me为使用等效构件法得到的计算长度系数[11]。
上文分别给出三阶、四阶、五阶、八阶阶梯柱的长度系数的数值计算结果，并且使用ANSYS 17.0对上述阶梯柱的组合进行了求解，并对本文递推公式得到的结果与ANSYS 17.0计算结果进行了对比，发现该递推公式与Levenberg-Marquardt算法相结合所求解出的长度系数具有很高的精度。
使用MATLAB2014b编制了n阶超越方程组的非线性数值解程序,该程序可以对9阶及9阶以下，长度和惯性矩任意组合的阶梯柱的临界应力和长度系数进行求解。
表5 八阶阶梯柱计算长度系数及比较
Tab.5 Calculated length factor and comparison of8 stepped column
β2/β3/β4/β5/β6/β7/β8μ(文献[11])μa(文献[9])Me(文献[11])数值法(2μ2)ANSYS17.0本文数值法与ANSYS 17.0方法相对误差1.2/1.2/1.2/1.2/1.2/1.2/1.22.451 72.451 72.449 72.413 12.413 101.3/1.3/1.3/1.3/1.3/1.3/1.32.739 22.739 22.731 42.680 02.679 90.003 71.6/1.6/1.6/1.6/1.6/1.6/1.63.909 93.909 93.812 23.782 43.782 20.005 31.3/1.3/1.3/1.3/1.2/1.2/1.22.705 32.705 32.699 22.648 72.648 701.6/1.6/1.3/1.3/1.3/1.2/1.23.248 13.248 13.229 03.172 33.158 70.430 61.8/1.8/1.6/1.3/1.3/1.2/1.23.863 33.863 33.817 63.759 73.740 80. 510 01.2/1.2/1.2/1.2/1.2/1.2/1.22.451 72.451 72.449 72.413 12.413 10
4 结论
本文使用伸缩臂的阶梯柱模型对n阶阶梯柱的稳定性微分方程组进行了推导，得到了n阶阶梯柱递推公式，根据阶梯柱模型的力学和结构特性，列写了补充方程；使用Levenberg-Marquardt数值算法求解n阶阶梯柱的超越方程组，与现行国家标准GB/T3811—2008和ANSYS 17.0所得结果进行对比，结果证明了n阶阶梯柱递推公式正确，数值解法相比其他的计算方法通用性更强，精度更高。此外，阶梯柱模型的长度系数具有一定的非线性，小范围内的插值不会产生太大的误差，但对于大截面的阶梯柱模型，若使用插值法计算，则临界力的误差较大。
[1] 陆念力，都亮，兰朋.变截面阶梯压杆精确失稳特征方程及其稳定计算实用方法——《塔式起重机设计规范》中阶梯柱计算长度的快捷计算及精度分析[J].建筑机械，2014(6):76-81.
LU Nianli， DU Liang， LAN Peng. Accurate Buckling Characteristic Equation of Stepped Column Bar and Its Stability Analysis with Practical Method：Shortcut Calculation and Accuracy Analysis with the Effective Length of Stepped Column Given in Specification for Tower Crane Design[J].Construction Machinery,2014(6):76-81.
[2] LI Wenjun, ZHAO Jiong, JIANG Zhen, et al. A Numerical Study of the Overall Stability of Flexible Giant Crane Booms[J]. Journal of Constructional Steel Research, 2015, 105: 12-27.
[3] TIMOSHENKO S P, GERE J M. Theory of Elastic Stability[M]. New York: Mc Graw-Hill, 1961: 49-63.
[4] 王欣，易怀军，赵日鑫，等.一种n阶变截面压杆稳定性计算方法的研究[J].中国机械工程，2014，25(13),1744-1747.
WANG Xin, YI Huaijun, ZHAO Rixin, et al. Research on Stability Analysis Method of n-order Variable Cross-Section Compression Bars[J]. China Mechanical Engineering,2014,25(13):1744-1747.
[5] 张宏生，陆念力.动臂式塔机变截面吊臂的整体稳定性分析[J].哈尔滨工业大学学报，2010，42(9):1394-1397.
ZHANG Hongsheng，LU Nianli.Global Stability Analysis of the Tapered Luffing Jib of Tower Crane[J].Journal of Harbin Institute of Technology, 2010，42(9):1394-1397.
[6] 刘庆潭.复杂细长压杆临界力计算的传递矩阵法[J].力学与实践，1997(6):51-54.
LIU Qingtan. Transfer Matrix Method for the Complex Compressed Bar [J]. Mechanics in Engineering,1997(6):51-54.
[7] 孙建鹏，李青宁，曹现雷.压弯杆弹性弯曲分析的精细传递矩阵法[J].北京工业大学学报，2012,38(4):536-539.
SUN Jianpeng，LI Qingning，CAO Xianlei. Precise Transfer Matrix Method for Elastic Bending Analysis of a Compressive-bending Bar[J]. Journal of Beijing University of Technology, 2012,38(4):536-539.
[8] 陆念力，兰朋，李良.二阶理论条件下的梁杆系统精确有限元方程及应用[J].哈尔滨建筑大学学报，1998,31(4):67-74.
LU Nianli, LAN Peng, LI Liang.Appilication of FEM of Beam Element with Theory of ⅡOrder[J].Journal of Harbin University of Civeil Engineering Architecture,1998，31(4)：67-74.
[9] 陆念力，兰朋，白桦.起重机箱形伸缩臂稳定性分析的精确理论解[J].哈尔滨建筑大学学报，2000,33(2):89-93.
LU Nianli, LAN Peng, BAI Hua. Precise Stability Analysis of Telescopic Boom[J]. Journal of Harbin University of Civil Enineering Architecture,2000，33(2):89-93.
[10] 陆念力，都亮.多级阶梯柱侧向刚度与轴压临界力的精确分析及其实用算式[J].工程力学,2015,32(8):217-222.
LU Nianli , DU Liang. Analysis of Lateral Stiffness and Critical Buckling Force of Stepped Columns and Their Practical Formulations[J].Engineering Mechanics ,2015,32(8):217-222.
[11] 中华人民共和国国家质量监督检验检疫总局，中国国家标准化管理委员会.GB/T 3811—2008起重机设计规范[S]. 北京：中国标准出版社, 2008.
General Administration of Quality Supervision, Inspection and Quarantine of the People’s Republic of China, Standardization Administration of the People’s Repubic of China. GB/T3811—2008 Design Rules for Cranes [S]. Beijing: China Standard Press,2008.
[12] 刘士明，陆念力，孟丽霞.牵绳非保向力作用下的起重臂稳定性分析[J].哈尔滨工业大学学报，2014,46(3):26-29.
LIU Shiming，LU Nianli，MENG Lixia.Stability Analysis of Telescopic Booms under Pull-rope Follower Force[J].Journal of Harbin Institute of Technology, 2014,46(3):26-29.
[13] 刘士明.工程起重机伸缩臂系统结构稳定性及复合运动动力学研究[D].哈尔滨：哈尔滨工业大学,2013.
LIU Shiming. Research on Structural Stability and Complex Motion Dynamic of Telescopic Boom System in Construction Crane[D].Harbin:Harbin Institute of Technology ,2013.
[14] BAGHERI S, NIKKAR A. Higher Order Explicit Solutions for Nonlinear Dynamic Model of Column Buckling Using Variational Approach and Variational Iteration Algorithm-Ⅱ[J]. Journal of Mechanical Science and Technology, 2014, 28(11): 4605-4611.
