刘宇新，倪长辉，张 洋，王雅妮

（1.中联重科股份有限公司建设机械关键技术国家重点实验室，湖南 长沙 410000；2.中联重科股份有限公司工程起重机分公司，湖南 长沙 410000）

细长梁柔性结构在各行各业中应用极为广泛，常见的如直升机螺旋桨、汽轮机低压末级叶片、卫星天线、起重机臂架、泵车折叠臂等。这些结构都具有相同的特点——几何非线性因素影响较大。

为提高细长柔性梁结构变形的计算精确度，从力学模型简化到计算分析方法已经进行了大量的研究。在模型简化方面，目前较为成熟的2 种简化模型分别为Euler 梁和Timoshenko 梁。常见的材料力学、结构力学等经典力学教程中，都是基于梁的小变形弹性假设进行处理的，该处理方法因未考虑几何非线性的影响，其计算结果误差较大。目前针对细长梁柔性结构变形的精准计算方法，国内外已有大量的研究。Saje［1］从基础理论出发，针对梁的变形分析提出了一种变分原理。关于梁的大变形分析方法，Dado 等［2］进行了系统总结，常见的有椭圆积分法［3］、打靶法［4］、有限差分法［5］及有限元法［6-7］。有限元法是目前应用最为广泛的一种方法，得益于计算机技术的发展。有限元法中最常用的2种方法分别为Euler描述法［8］和Lagrangian描述法［9］。

随着人工智能（artificial intelligence，AI）技术的发展，机械行业对智能控制等AI 技术需求相应增加，对长细梁柔性结构精准控制时，需要精确的变形响应计算方法，并且响应时间要求在毫秒级。传统的有限元及非线性迭代计算方法虽然精度能够满足响应要求，但计算时间不能满足毫秒级要求，因此急需一种满足一定精度要求的快速变形计算方法。本文以有限元及梁偏微分方程为基础，推导一种快速计算梁变形的计算方法。

1 梁的基本方程

梁理论对于挠曲线的微分方程为［10］

式中：ρ为梁挠曲线的曲率半径；θ为梁挠曲线x处的截面转角；s为挠曲线弧长；M（x）为梁的弯矩方程；E为材料弹性模量；I（x）为梁横截面惯性矩。

小变形时，θ很小，可以近似认为：tanθ=sinθ=θ，并且忽略了梁的轴向变形的影响。常见的等截面悬臂梁端部受力结构如图1所示。

图1 悬臂梁载荷关系Fig.1 Schematic of cantilever beam loads

梁的微分控制方程为［11］

解该方程通常需要采用数值积分的方法进行求解，需要不断进行迭代求解，相对计算时间较长。同时对于变截面连续梁，其微分形式比式（2）更为复杂，其求解难度更大。

2 工程算法探索

悬臂梁作为工程机械行业中常见的机械结构形式，在工作过程中的变形状态是从设计到使用阶段中影响其性能的一个重要参数指标。在设计阶段，对变形的准确计算能够有效地指导设计，从而设计出性能更优的悬臂梁结构；在使用阶段，对变形的准确计算便于对机械进行精准智能控制，是自动控制技术必不可少的基础。

对于悬臂结构目前已知的几种典型变形计算方法包括：放大系数法、叠加法、等效惯性矩法、微分方程数值求解等解析法和有限元法（小变形、大变形）。

以上6 种常用悬臂梁变形计算方法基本囊括了目前处理悬臂梁结构的计算方法，并且这些计算方法的共同输入参数有：梁的长度L，梁的截面惯性矩I，梁的弹性模量E，梁的外载荷P。其中梁的长度及惯性矩还需要根据每一段梁的不同截面分别给出。最为重要的一点是以上各分析方法都是建立在简化力学模型的基础上完成的，在计算时未考虑制造误差、装配间隙等因素对最终变形计算结果的影响。因为目前的计算需求通常在设计阶段，一般采用安全系数的方式覆盖这些因素的影响，只需要保证在安全系数覆盖的前提下满足相关设计规范的要求即可。

随着精准控制要求越来高，单纯的只是保证结构安全已经不足以满足需求。为了实现实时的精准控制，需要对悬臂结构工作后变形量进行精确判定。为了做到变形量的实时动态输出，其计算响应时间通常要求小于500 ms，同时计算精度也需要满足需求，这对结构变形计算提出了新的要求。传统计算变形的方法已经不能同时满足精度及计算响应时间要求。

本文以梁微分方程为基础，采用泰勒展开的方式，通过结合实测数据进行数值拟合的形式建立梁的变形解析表达式。泰勒展开数学描述为

若函数f(x)在包含x0的某个区间（a，b）内具有（n+1）阶导数，则对于任意x∈(a，b)，有

结合梁变形后曲线及式（2），在（0，L）内存在多阶导数，对于悬臂梁结构的变形方程y(x)，x∈[0，L]，令x0=0，显然y(x0)=0。则对于悬臂梁结构其变形方程可以根据式（3）结果简写为（忽略高阶项，这里只保留3阶，且C0=0，C1=0）：

对于式（4）中的系数项，由式（2）很难直接求出其具体数值，这里2个未知数C2和C3，需要至少2个已知变量才能求出式（4）的具体表达式。

由于以上推导都是基于力学简化的梁模型获得的，实际梁结构中还存在制造误差、装配间隙等因素的影响，因此，要获得较为准确的式（4）表达式，需要通过与试验数据相结合的方式确定最终的未知参数C2和C3。对于实际产品，大多数产品出厂前会进行强度、刚度（变形）的验证试验，此时可以获得产品的真实变形量，将这些变形量作为已知条件代入式（4）则可以求出具体的表达式。

3 算例

本文以汽车起重机的臂架结构为例进行数值验证，某汽车起重机主臂结构如图2所示。由于变幅油缸的轴向刚度很大，变幅油缸与伸缩主臂铰接点以下及起重机工作台形成一个稳定的三角形区域，臂架铰点到变幅点之间基本不会产生变形，因此可以将主臂实际结构等效为根部固定的悬臂梁。

图2 汽车起重机主臂结构Fig.2 Schematic of boom

由图2可知：以吊臂铰点为原点、悬臂梁轴向为x轴、悬臂梁横向为y轴，建立坐标系。将吊载沿坐标系分解为横向力T和轴向力N。若采用传统解析法中计算效率最高的等效惯性矩法求解，则端部挠度计算公式描述为

式中：f1=为横向力产生的吊臂端部挠度；Cf为放大系数［12］；Ncr为临界力；E为弹性模量；I为等效惯性矩。

在求解等效惯性矩与临界力时，计算公式复杂，需输入每节臂的臂长、惯性矩及臂节间搭接长度等参数。

使用本文实测数据结合数值拟合方法，只需要代入不同工况下的臂头挠度实测数据，即可反求得等效惯性矩I、临界力Ncr等参数，从而得到吊臂端部挠度的具体表达式，在此基础上，采用数值拟合方法，即可得到整个吊臂的挠度变形曲线。

以500 t 全地面起重机主臂全伸工况为数值算例，展开验证分析，主臂参数信息见表1。为了描述方便，对主臂的面积和惯性矩数值进行归一化处理，设第一节臂的面积和惯性矩为1。

表1 臂架参数Tab.1 Boom parameters

如图3所示，基于ANSYS 按表1实际尺寸沿吊臂轴线方向建立吊臂的有限元仿真计算模型。约束吊臂后铰点，将吊载沿坐标系分解，打开结构的大变形效应。在主臂全伸工况下，选取了2 种不同仰角，每个仰角选取了3 种不同吊重。6 种工况组合见表2。

图3 汽车起重机主臂结构有限元模型Fig.3 Finite element model of boom

表2 工况列表Tab.2 Tableofload cases

在全伸臂长L=85 m，设EI=a，Ncr=b，则式（5）可以转化为

代入2个不同工况，即可求解参数a、b，从而得到吊臂端部挠度的表达式。因此，将工况1和工况2载荷和挠度数据，代入式（6），得到a=9.67×109，b=1.76×106。根据式（6）及L=85 m、a=9.67×109、b=1.76×106，可以求得工况3～6的吊臂端部数据，并将该数据与ANSYS仿真数据对比，误差见表3。由表3可知，所有工况的吊臂端部挠度绝对误差不超过0.25 m。

表3 吊臂端部挠度数据误差对比表Tab.3 Comparison table of displacement error onboom ends

基于吊臂端部挠度数据及式（4）悬臂梁挠曲线方程，将吊臂等效为等截面梁进行挠曲线数值拟合，数值拟合曲线方程的解析式为

式中：a=9.67×109；b=1.76×106；L=85。

工况1～6工况的数值拟合曲线与仿真变形曲线对比如图4所示。由图4可知，6个工况下2条曲线的挠度最大差值不超过0.5 m。对于500 t全地面起重机，在主臂全伸工况下，起重量不大于20 t。因此，本文的工程算法即实测数据结合数值拟合方法具有很强的工程适用性，且误差在工程允许范围内。

图4 挠度曲线对比Fig.4 Comparison diagram of displacement curves

4 结语

本文采用实测数据与数值拟合相结合的方式，能够快速、准确地实时求出工程机械中梁结构的变形数学表达式，这与传统的基于力学简化模型的计算方法相比，避免了输入参数多的情况，在计算效率与计算精度无法同时满足需求的情况，为实时智能控制提供了保障。