基于格点型有限体积法的活塞高周疲劳分析

龚京风，徐宗著，宣领宽

(1.武汉科技大学汽车与交通工程学院，湖北 武汉，430065；2.中国舰船研究设计中心船舶振动噪声重点实验室，湖北 武汉，430064)

活塞是柴油机的核心部件，其长期处于高温、高压的工作环境且承受应力水平较高的周期性载荷，活塞的抗疲劳性能直接影响到柴油机的工作可靠性。活塞在工作过程中的疲劳损伤包括高周疲劳损伤和低周疲劳损伤，前者主要是指柴油机因高频变化的温度、燃气爆发压力及往复惯性力而导致的损伤，后者主要是指柴油机在启动、停机等工况变化中因活塞温度大范围变化引起的塑性变形和蠕变等损伤[1]。由于柴油机长时间处于工作循环中，而启停工况下的运行时间较短，因此本文主要对活塞高周疲劳进行分析。

在计算力学领域，除了FEM之外，有限体积法(FVM)也得到了广泛关注，尤其是格点型有限体积法(CV-FVM)。Xia等[8]将CV-FVM应用于二维和三维悬臂梁的结构动力学计算，所得到的解与解析解一致。陈浩[9]基于确定性结构动响应的控制方程推导出随机结构动响应的控制方程，并使用CV-FVM对控制方程进行离散，给出了随机结构动响应控制方程的求解过程。井丽龙[10]通过直梁静力弯曲分析、固有频率分析以及集中载荷动力响应分析，验证了CV-FVM在梁、板、壳这类结构的静力学和动力学分析中的准确性。笔者等[11-13]使用CV-FVM对热弹性问题进行了研究：在文献[11]中利用双线性四边形单元和三角形单元的混合网格分析功能梯度材料的热弹性问题，发现相比于梯度单元和高阶单元，CV-FVM可以减少计算量和内存消耗；在文献[12]中研究了各向异性材料的热弹性问题，拓展了CV-FVM的应用范围；在文献[13]中研究了考虑温度热效应的功能梯度材料的热弹性问题，进一步验证了CV-FVM的计算精度和适用性。宣领宽[14]在研究结构声耦合问题时采用三角形单元、四边形单元、四面体单元、三棱柱单元和六面体单元划分计算域，基于CV-FVM离散各向异性材料的弹性力学控制方程，分别使用显式和隐式算法求解离散方程，并将该方法得到的结果与商用软件计算结果进行了对比。

1 高周疲劳数值模型的建立

1.1 CV-FVM控制体简介

以图1中节点N1为例对CV-FVM控制体进行说明。C1、C2分别为N1的两个相邻单元1和2的中心；F1、F2、F3、F4、F5分别为N1相邻网格面单元的中心，L1、L2、L3、L4分别为N1相邻网格单元边线的中心。依次连接这些中心点，即为N1的控制体，其中，四边形C2F4L3F5、C2F4L4F3、C2F5L2F3为N1在单元2内的控制体的边界面，四边形C1F1L1F2、C1F5L3F1、C1F2L2F5为N1在单元1内的控制体的边界面。

图1 CV-FVM控制体示意图

1.2 弹性方程的离散

设控制体的体积为V、边界面积为S，根据力平衡条件建立弹性方程：

(1)

式中：u、v、w分别为坐标轴x、y、z方向的位移分量；nx、ny、nz分别为控制体边界面的法向矢量在x、y、z方向的分量；fx、fy、fz分别为体积力在x、y、z方向的分量；ρ为密度；σ表示应力。

由文献[15]中的本构方程和几何方程得

(2)

将式(2)代入式(1)中，考虑给定面力边界条件的影响，得

(3)

(4)

(5)

式(3～5)中：pxB、pyB、pzB分别为x、y、z方向的面力；SN为边界面的面积。当给定位移边界uiB(i=x,y,z)时，边界上的节点位移满足ui=uiB，可采用文献[16]中的置大数法处理给定位移边界。

以式(3)为例，对弹性问题的控制方程进行离散，将该式右端第一项和第二项离散为

(6)

式中：nc为当前节点的相邻单元总数；ncni为当前节点第i个相邻单元上的节点总数；Si为当前节点在相邻单元i内的控制体边界面积；Nij为当前节点第i个相邻单元上第j个节点处的形函数。关于不同单元类型的形函数导数的积分计算方法见文献[14]。

将式(3)右端第三项离散为

(7)

式中：Vi为当前节点第i个相邻单元的体积；fxi为当前节点第i个相邻单元在x方向的体积力。

将式(3)右端第四项离散为

(8)

式中：nN为当前节点相邻边界网格面的个数；ABi为当前节点在第i个相邻边界网格面内的控制体边界面积。

对式(3)左端项采用中心差分法离散为

(9)

本文主要解决静态弹性问题且不考虑体积力，无需考虑时间项和体积源项，因此式(3)的离散方程可简化为

(10)

式(4)和式(5)的离散方程可用相同方式得到：

(11)

(12)

将计算域划分为具有n个节点的网格模型，对每个节点联立方程式(10)～(12)，得到一个求解弹性问题的线性方程组，表示为

AX=B

(13)

其中

A=

式中：列向量B由方程中已知边界和体积力求得，其中下标un、vn、wn表示各位移方向上不同节点对应的已知系数；系数矩阵A中各节点上未知变量的系数由形函数导数积分以及部分给定位移边界求得，其中下标un1、vn1、wn1表示各位移方向不同节点上的第一个方程对应的未知量系数，un2、vn2、wn2、un3、vn3、wn3定义类似。

由几何方程以及形函数可求得各节点上的应变分量：

(14)

再由式(2)计算得到计算域内各节点的应力分量，同时可由式(15)求得von Mises等效应力

(15)

1.3 高周疲劳寿命计算方法

在疲劳寿命分析中，对于恒幅应力循环载荷，描述载荷谱的参数主要有最大应力σmax和最小应力σmin，同时可以根据σmax和σmin求得应力幅σa、平均应力σm和应力比R，当R=-1时，称为对称循环载荷。

高周疲劳寿命分析时，通常采用S-N曲线来表征应力幅与疲劳寿命之间的关系，恒幅应力循环载荷作用下，采用应力比和应力幅描述循环应力水平，当已知应力比时，应力幅即为主要评价参量。描述材料疲劳性能的基本S-N曲线大多是根据试验得到的S-N数据通过拟合得到曲线表达式。本文使用分段线性插值，即每两个数据点使用线性函数拟合，得到S-N曲线表达式(16)：

σa=f(N)

(16)

式中：N为材料疲劳寿命。

材料应力和疲劳寿命数据都是利用光滑小尺寸试件在对称循环载荷作用下得到的，对于大尺寸构件，需要引入一个小于1的尺寸修正因子；类似地，表面粗糙度的影响也可以引入一个小于1的修正因子。将这两种影响因素综合考虑，定义为疲劳强度因子kf，代入到式(16)得

σa=kff(N)

(17)

文献[17]指出，在相同的应力幅下，随着平均应力的上升，构件的疲劳寿命将缩短，因此在进行疲劳寿命分析时需要考虑平均应力的影响。可采用以下两种方法对应力幅进行修正。

(1)使用Goodman理论修正：

(18)

(2)使用Gerber理论修正：

(19)

相比较来看，平均应力增大到一定程度后，Goodman理论修正的应力幅值要高于Gerber理论修正的应力幅值，即Goodman理论比Gerber理论在疲劳寿命预测上更加保守。在工程实际中往往采用Goodman理论对应力幅值进行修正。

对于静态弹性问题，给定最大外载荷，将计算得到的σvon_Mises作为载荷谱最大应力，再由给定的应力比计算得到应力幅，即可对构件的疲劳寿命进行预测。

2 数值模型的验证

在板左侧分别施加90、100、110 MPa的法向压力载荷。板右侧和上侧设置为对称边界，下侧为自由边界，并指定应力比R=0.1，kf=1，使用Goodman修正理论进行应力幅修正，验证模型的S-N曲线见图3。

图2 验证模型

图3 验证模型的S -N曲线

分别使用FEM和CV-FVM方法计算该模型的应力，并将最大等效应力计算结果与文献[18]中的研究结果进行对比，见表1。

表1 最大等效应力计算结果

从表1可以看出，在3种压力载荷作用下，CV-FVM法与其他两种方法的计算结果误差在1%以内，验证了CV-FVM方法的正确性。分别使用FEM和CV-FVM方法计算在110 MPa法向压力载荷作用下矩形板的疲劳寿命，结果如图4所示。

(a) FEM计算结果

(b) CV-FVM计算结果

使用FEM和CV-FVM方法计算得到的最小疲劳寿命分别为221 000和226 758循环周次，两者数值极为接近，疲劳寿命云图分布也呈现出相同规律，这充分验证了本文提出的高周疲劳数值模型的正确性。

3 活塞高周疲劳寿命分析

3.1 活塞计算模型及边界条件

活塞材料为ZL109铝合金，密度为2770 kg/m3，杨氏模量为69 GPa，泊松比为0.33。忽略倒角等几何特征，采用四面体单元划分计算域。活塞计算模型见图5，该模型为完整活塞的1/4模型。

图5 活塞计算模型

研究活塞疲劳问题时，活塞受到的燃气压力载荷可以近似为正弦脉冲循环载荷[19]。本文在给定最大燃气压力Pg和应力比的基础上，基于静态应力结果模拟计算正弦脉冲循环载荷作用下的活塞疲劳寿命。在活塞顶部、火力岸及第一道环上表面施加157.3 MPa的燃气压力，第一环槽内侧及下表面按燃气压力的75%施加，第一环岸、第二环槽上表面及下表面按燃气压力的25%施加，第二环槽内侧按燃气压力的20%施加[20]。同时还考虑此时的支反力FN作用在活塞销孔上表面的表面压力PN，支反力与表面压力的关系满足下式：

(20)

式中：l为销座长度；r*为销孔半径。

所研究的活塞中，一侧销孔的支反力FN=177 kN，l=59.7 mm，r*=65 mm，可得PN=12.83 MPa。为约束活塞的刚体位移，在两个剖面分别施加对称边界条件，销孔上表面约束活塞轴线方向位移。

活塞疲劳计算中给定应力比R=0.01，极限强度σu=239.3 MPa，S -N曲线上部分数据见表2。

表2 活塞S -N曲线上部分数据

活塞尺寸与材料疲劳性能试件的尺寸存在较大差异，需考虑尺寸效应，对于圆截面构件，可由式(21)计算疲劳强度因子[17]。活塞在加工过程进行了抛光处理，表面粗糙度修正因子为1，即可忽略表面粗糙度对活塞疲劳寿命的影响[21]。

(21)

式中：d为活塞直径，本文中d=0.0835 m。计算得到kf=0.77。

3.2 活塞高周疲劳计算结果

按前述边界条件，分别使用FEM和CV-FVM计算活塞应力场，得到活塞的总位移云图(见图6)。从图6可以看出，两种方法计算得到的云图在活塞位移分布规律上高度吻合。活塞最大位移主要集中在头部区域，活塞销孔位置由于y方向的位移被约束，最小位移出现于此处。

(a) FEM计算结果

(b) CV-FVM计算结果

图7为两种方法计算得到的等效应力云图，其中FEM计算得到最大等效应力为86.1 MPa，CV-FVM计算得到的最大等效应力为86.5 MPa，两种方法的计算结果均在同一范围内变化，云图分布规律也一致，最大等效应力主要集中在活塞销孔上表面内侧区域。

(a)FEM计算结果

(b) CV-FVM计算结果

在销孔上侧边沿取特征点(见图8)上的等效应力值，比较两种方法在数值上的精度，见图9。从图9中可以看出，CV-FVM和FEM计算得到的等效应力数值吻合良好，验证了使用CV-FVM计算活塞应力场的准确性。

图8 特征点示意图

图9 特征点上的等效应力

在上述计算基础上，使用Goodman和Gerber平均应力修正理论计算活塞高周疲劳寿命，结果分别为8.88×1014循环周次和3.69×1015循环周次。为保证构件的使用安全，一般对疲劳寿命进行保守估计，因此这里使用Goodman修正方法的计算结果。图10分别为FEM和CV-FVM计算得到的疲劳寿命云图。FEM所得最小疲劳寿命为1.01×1015循环周次，CV-FVM所得最小疲劳寿命为8.88×1014循环周次，二者基本一致。在机械载荷的作用下，活塞最小疲劳寿命出现在销孔上表面内侧区域，与最大等效应力位于同一区域，即活塞等效应力与疲劳寿命成反比，考察机械载荷作用下的活塞疲劳破坏，要重点关注其销孔处的应力情况。

(a)FEM计算结果

(b) CV-FVM计算结果

4 结语

使用本文数值模型预测活塞高周疲劳寿命，其结果与有限元法计算结果误差很小。在机械载荷作用下，活塞等效应力与疲劳寿命成反比，最小疲劳寿命出现在活塞销孔上表面内侧区域，与最大等效应力出现的区域一致。