抽油机是开采石油的一种机器设备，俗称“磕头机”，在石油开采行业被广泛应用，是有杆抽油系统最主要的举升设备，根据其是否有游梁，可分为游梁式抽油机和无游梁式抽油机[1]。一台抽油机设备通常有几十吨甚至上百吨的重量，其支座结构占据其较大比重，是抽油机的重要受载结构，为整个设备运转以及各零件组合运动提供支撑作用。因此，对抽油机支座结构进行优化设计至关重要。
现代设计中，轻量化设计已经成为优化设计的主要工作，可在实际生产与工程作业中降本增效。马天宇等[2]在渭河大桥设计参数优化问题中，提出加劲钢桁架桥梁设计参数化优化方法，借助有限元软件Midas/Civil建立模型，使主跨用钢量降低2.14%。张明康等[3]在中低速磁悬浮列车悬浮架纵梁结构优化问题中，分别基于多项式响应面法、Kriging法构造目标参数关于设计变量的代理模型，利用多目标优化算法进行寻优，最终得到满足强度可靠性要求下，减重16.9%的纵梁结构。方志凌等[4]在汽车前机舱一体化压铸设计问题中，通过固体各向同性材料惩罚模型法(solid isotropic material with penalization，SIMP)对前舱进行拓扑优化，再经过仿真分析，最终成功使产品减重13.9%。可见，轻量化设计可以大幅减少材料用量，节约成本。
本文以最新设计的12-8CYJ-00型抽油机为研究对象，在已有的设计基础上对其支座部分进行轻量化设计。广泛使用的结构优化方法有结构拓扑优化、结构参数优化以及结构形状优化。结构拓扑优化方法会使现有的几何拓扑结构发生较大变化，将增加切割、冲孔等减材制造工艺，反而增加生产成本，而结构形状优化在早期设计中已经被充分考虑，因此本次轻量化设计采用结构参数优化的方式实现。
本文利用ABAQUS有限元软件平台搭建一套结构优化框架，并完成轻量化设计工作。通过有限元软件建立参数化模型，再基于模型的有限元计算结果与结构设计参数之间的关系，采用灵敏度分析对支座结构进行优化，要求优化后的抽油机支座在任意工况下满足强度条件的同时，结构整体质量也更轻量化。
1 强度分析及参数化建模
1.1 工程背景
12-8CYJ-00型抽油机结构如图1所示。其中，工作臂采用四连杆结构，由电动机驱动，工作臂前端部负载12 t，尾部为了平衡负载，挂有12 t的配重块。结构总重量（包含负载与配重块）为58 t，其中支座部分（包含前架、后架、加强板）的重量为16.6 t。支座结构的材料为Q235R碳素结构钢，材料的弹性模量E、泊松比ν、密度ρ、屈服强度σs以及许用应力[σ]如表1所示，支座结构整体由3种型号的矩形空心钢构成主要支座部分，截面如图2所示。在后架处设计有20 mm厚的加强板。设备用于两种工况，具有不同的运转速度，工况1和工况2下设备运转一圈的运动周期分别为30、15 min。要求对支座结构进行减重设计，同时保证其强度要求，采用Mises应力准则作为强度准则。
<G:\武汉工程大学\2025\第4期\宋永军-1.tif>[四连杆组件][负载悬挂][驱动安装架][前架][混凝土基座][加强板][配重块][动力组件][后架][工字钢基座]
图1 抽油机结构
Fig. 1 Structure of oil pumping rig
表1 Q235R材料参数
Table 1 Material properties of Q235R
[E / MPa ν ρ / （kg·m-3） σs / MPa [σ] / MPa 210 000 0.3 7.9×103 235 200 ]
<G:\武汉工程大学\2025\第4期\宋永军-2.tif>[10][10][10 ][unit： mm][200][200][300][100][200][200][Ⅰ Ⅱ Ⅲ]
图2 各型号矩形空心钢截面
Fig. 2 Rectangular hollow steel sections of each type
1.2 支座荷载提取与边界条件建立
1.2.1 支座荷载提取 在设备运转过程中，连杆组件与支座上的反力保持动力学平衡，使用ABAQUS软件平台完成有限元计算，得到连杆机构对支座的反力。将模型分为连杆组件和支座两大部分，先对连杆组件进行计算，其关节转动副采用ABAQUS软件的铰接连接器单元替代，底部与尾部设定参考点耦合，对两个位置施加对应的位移约束，如图3所示。
<G:\武汉工程大学\2025\第4期\宋永军-3.tif>[12 t][12 t][Y][Z][O]
图3 连杆组件的边界条件
Fig. 3 Boundary conditions for the connecting rod assembly
对上述两种工况下的一个运动周期进行仿真计算，从计算结果中对连杆组件底部的反力进行提取，将两种工况的运动周期归一化，得到两种工况下前架反力在一个运动周期的变化，如图4所示。提取两种工况下连杆组件尾部反力作用位置的Z坐标及两种工况下的后架反力，在归一化后的运动周期内的变化如图5所示。
<G:\武汉工程大学\2025\第4期\宋永军-4.tif>[0.0 0.2 0.4 0.6 0.8 1.0
运动周期][600
500
400
300
200
100
0][反力 / kN][Z方向反力
Y方向反力][（b）][（a）][Z方向反力
Y方向反力][0.0 0.2 0.4 0.6 0.8 1.0
运动周期][800
600
400
200
0][反力 / kN]
图4 前架反力运动周期变化：(a) 工况1，(b) 工况2
Fig. 4 Periodic change of the reaction force of the front
support：(a) case1，(b) case2
<G:\武汉工程大学\2025\第4期\宋永军-5.tif>
图5 后架反力及连杆尾部Z坐标随运动周期变化
Fig. 5 Periodic change of the reaction force of the rear
support and the Z coordinate of the tail
1.2.2 支座边界条件建立 在ABAQUS软件中建立模型，其主要结构由梁单元与壳单元组成。根据得到的各工况下的连杆组件反力变动情况，对支座模型施加如图6所示的边界条件，前架载荷幅值取决于连杆组件底部对前架顶部的反力，而后架上的载荷幅值以及作用位置都取决于连杆组件尾端对后架的反力。支座底部采用固定的位移边界条件。
<G:\武汉工程大学\2025\第4期\宋永军-7.tif>[O][Z][Y][分布力
集中力
固定
]
图6 支座边界条件
Fig. 6 Boundary conditions for the support
对于后架上的载荷变动情况，需要通过ABAQUS软件提供的Dload子程序[5]定义。用Fortran编程语言编写Dload子程序需要有相应变量的解析式，由于后架载荷既有幅值上的变动，又有位置上的变动，二者都是时间t的函数，因此对后架上的载荷变动采用分段线性插值的方式模拟，其结果如图7所示，得到的载荷变动表达式如式(1-3)所示，其中P(t)是反力Z坐标随运动时间变动的函数，F1(t)、F2(t)分别是工况1和工况2下载荷幅值随时间变动的函数。
<G:\武汉工程大学\2025\第4期\宋永军-8.tif>
图7 后架荷载与位移多段线性插值
Fig. 7 Multi-segment linear interpolation of rear support load and displacement
[P(t)=19 655t+6 505.5,-14 601.4(t-0.11)+8 667.5,-6 692.5(t-0.40)+4 396.6,795.9(t-0.6)+2 041.5,16 527.2(t-0.8)+3 200.5,0≤t≤0.110.11<t≤0.400.40<t≤0.600.60<t≤0.800.80<t≤1] (1)
[F1(t)=-735 591t-45 914.5,-1 800 096(t-0.05)-82 496,1 198 572.2(t-0.11)-194 075,432 774(t-0.2)-86 203,136 957(t-0.3)-42 926,13 298.6(t-0.4)-29 230.4,-162 865.7(t-0.86)-23 113,0≤t≤0.050.05<t≤0.110.11<t≤0.200.20<t≤0.300.30<t≤0.400.40<t≤0.860.86<t≤1]
(2)
[F2(t)=-1 090 880t-45 286.7,-3 445 266.7(t-0.05)-99 830,2 497 844(t-0.11)-306 546,413 950(t-0.2)-81 740,127 970(t-0.3)-40 345,16 080.4(t-0.4)-27 548,-3 090 827(t-0.86)-20 151.1,0≤t≤0.050.05<t≤0.110.11<t≤0.200.20<t≤0.300.30<t≤0.400.40<t≤0.860.86<t≤1]
(3)
1.3 支座参数化模型建立
两种工况在ABAQUS中进行有限元分析，得到Mises应力云图，如图8所示，工况1下最大Mises应力为 68.558 MPa，位于后架处，工况2下最大Mises应力为181.461 MPa，同样位于后架。
<G:\武汉工程大学\2025\第4期\宋永军-9-1.tif><G:\武汉工程大学\2025\第4期\宋永军-9-2.tif>[Mises应力 / MPa][68.558
62.845
57.132
51.418
45.705
39.992
34.279
28.566
22.853
17.139
11.426
5.713
0.000][max][min][工况1][Mises应力 / MPa][max][min][181.461
166.339
151.218
136.096
120.974
105.852
90.731
75.609
60.487
45.365
30.244
15.122
0.000][工况2]
图8 优化前支座Mises应力云图
Fig. 8 Optimized the Mises stress contour diagram of the front support
根据两种工况下边界条件以及计算结果的对比可知，工况2对结构强度的要求高于工况1。因此选择工况2的边界条件作为参数优化模型的边界条件。
支座部分设计以板与桁架结构为主，其中加强板的设计变量主要由加强板厚度d1决定。支座主结构由3种型号的矩形空心梁通过焊接形成，主要设计参数有矩形截面的宽度、高度与截面壁厚，前架部分使用Ⅰ型截面梁，后架部分使用Ⅱ型截面梁与Ⅲ型截面梁，3种型号截面梁的截面及加强板的结构参数与取值范围见表2。
表2 结构参数表
Table 2 Structural parameters
[截面参数 设计变量 取值范围 / mm 加强板厚度 d1 8~20 Ⅰ型截面宽度 b1 60~100 Ⅰ型截面高度 h1 100~200 Ⅱ型截面宽度 b2 150~200 Ⅱ型截面高度 h2 150~200 Ⅲ型截面宽度 b3 150~200 Ⅲ型截面宽度 h3 250~300 截面壁厚 d2 6~10 ]
在ABAQUS软件中，采用梁单元建立支座部分桁架结构，采用壳单元建立加强板。二者使用绑定约束关系，该模型包含8个设计变量，即8个截面参数，这些变量通过在ABAQUS的inp文件中使用*Parameter关键字实现参数化[6]，成为参数化设计变量。
2 基于灵敏度分析的参数优化
2.1 优化目标函数与约束条件
结构的整体质量m取决于8个设计变量，因此采用单目标优化，优化目标函数如下：
[minm=2ρd2i=13Ai(bi+hi)+ρA4d1] (4)
式中：Ai为结构固定参数，i取1~4，A1、A2、A3分别为Ⅰ、Ⅱ、Ⅲ型截面梁使用长度，A4为加强板的面积，A1、A2、A3分别为73 015、25 900和94 948 mm，A4为24 749 891 mm2。
式(4)两边同时除以密度ρ，则新的优化目标函数为体积V，如式(5)所示。
[minV=2d2i=13Ai(bi+hi)+A4d1] (5)
当边界条件确定，通过有限元计算得到的最大Mises应力[σmax]可以看作是上述设计变量的函数[f(X)]，其中X为支座结构的设计变量所组成的向量，[X=[d1,b1,h1,b2,h3,b3,h4,b4]]，其形式如式(6)所示：
[σmax=f(X)] (6)
[σmax<[σ]=200 MPa] (7)
2.2 Python-ABAQUS优化求解框架
灵敏度分析是分析输出响应对设计变量变化的敏感性，灵敏度的数值可反映出各设计变量对输出响应变化的贡献率[7]。灵敏度分析在工程设计和概率安全评价中有着广泛的应用，其目的是探究系统或模型输出不确定性向输入不确定性的分配问题[8]。
灵敏度分析在优化设计问题中被广泛运用。在结构动态响应优化[9]、几何优化、参数优化以及系统优化[10-12]等工程问题中效果显著。
本次研究搭建了一个Python-ABAQUS优化求解框架[13]，适用于板壳、桁架结构减重设计。通过Python实现优化算法，调动ABAQUS隐式求解器，解算inp文件，通过有限元模型计算最大Mises应力来获取输出响应[σmax]。对函数[f(X)]的每一个变量求偏导数，即可获得8个设计变量对结构强度的灵敏度，如式(8)所示。
[Sj=?f(X)?Xj] (8)
式中：[Xj]为向量X的第j个分量，[Sj]是第j个分量的灵敏度，j取1~8。
整个优化流程是先求出结构强度对各个设计变量的灵敏度，用于衡量设计变量对结构强度的贡献率，对于贡献率低的设计变量，优先进行优化，再取一定的步长h，让该设计变量按此步长进行搜索，搜索方向为该设计变量负增长方向，步长h取值越小，越接近最优解。每优化一次，都会计算当前结构中最大Mises应力对设计参数的灵敏度，用来重新确定优化参数的优先级。灵敏度的数值求解方法是对各参数取一个微小扰动值δ，计算结构最大Mises应力变化与扰动值δ之间的比值。微小扰动值的选取依据为：当扰动值取2δ时的灵敏度与取值δ时的灵敏度之间相对误差小于2%即可认为扰动取值选取足够小。优化流程如图9所示。
2.3 参数优化与结果
将抽油机支座结构的初始设计变量输入优化程序中，Python主程序自动将参数代入ABAQUS的inp文件计算参数中，通过os. system()函数调用ABAQUS的standard.exe求解器与Dload子程序进行模型计算，主程序再得出结构最大Mises应力并计算对相关参数的灵敏度，通过灵敏度来确定参数优化的优先级。步长h的选取要求是当某一设计变量减小h，结构最大Mises应力的变化不超过10%，经过调试，取合适步长h = 5对优化参数进行搜索，反复进行迭代计算。每次计算完成，都要通过强度约束条件与参数边界约束条件判断，若强度达到结构许用应力200 MPa，则程序停止迭代并输出优化结果。若某一参数达到其边界值，则在下一次迭代中将其剔除，对其余变量进行优化，当全部变量都已达到边界约束时，程序停止运行并输出优化结果。
经过优化后，目标函数的结果为1.83 m3，为优化前的87.56%，即优化使得结构质量减少12.44%，总重量由原来的16.511 t减少到14.457 t，减重2.054 t。目标函数与最大Mises应力随优化迭代次数的变化见图10，优化后的参数如表3所示。
<G:\武汉工程大学\2025\第4期\宋永军-11.tif>
图10 目标函数与最大Mises应力随迭代次数的变化
Fig. 10 Changes in the objective function and the maximum Mises stress with the number of iterations
表3 优化前后参数对比
Table 3 Parameters before and after optimization mm
[设计参数 优化前 优化后 b1 100 60 h1 200 125 b2 200 150 h2 200 150 b3 200 180 h3 300 300 t1 20 20 t2 10 10 ]
对优化后的结构进行强度校核，结果如图11所示。在ABAQUS软件中计算得到其最大Mises应力为 199.629 MPa，位于支座结构的后架处，相较于优化前增长10.01 %，已经达到许用应力极限。
<G:\武汉工程大学\2025\第4期\宋永军-12.tif>[Mises应力 / MPa][max][min][199.629
182.994
166.358
149.722
133.086
116.450
99.815
83.179
66.543
49.907
33.272
16.636
0.000]
图11 优化后的Mises应力云图
Fig. 11 Optimized Mises stress contour
3 结 论
(1)优化前后结构参数对比得出，8个结构参数中，有3个未被优化，其余5个参量均已被优化。被优化的结构参数分别是Ⅰ型截面参数b1、h1，Ⅱ型截面参数b2、h2及Ⅲ型截面参数b3；未被优化的结构参数分别是矩形空心梁截面的厚度t2、加强板厚度t1及Ⅲ型截面参数h3。结果表明未被优化的参数，对最大Mises应力非常敏感。
(2)通过计算结构强度对设计参数的灵敏度，得出了在即时构型下各设计参数对结构强度的贡献率，对贡献率低的参数进行优化，并搭建一套Python-ABAQUS的优化框架，该框架便捷高效，适用于板壳、桁架结构的参数优化。经过优化的支座结构相较于优化前减重12.44%，减重2.054 t，优化效果显著。
(3)在对支座载荷的分析工作中，支座结构的前架与后架的距离对载荷有显著影响，若对其进行调控，有望减小载荷变动的峰值，改善抽油机的负载强度从而进一步对整个抽油机进行优化。