实际工程设计都必须考虑缺口特征对构件疲劳损伤的影响[1]，在相同的应力集中系数下，不同缺口特征导致不同的应力分布，即缺口效应[2]．缺口的存在不仅引起应力集中，还会导致缺口根部出现多轴应力场[3]，诱导更为复杂的应力应变响应，疲劳裂纹更易在这些区域萌生和扩展[4]．随着对应力集中研究的深入，由最初考虑名义应力发展到考虑局部应力应变，并最终认识到应力集中区域的应力梯度[5]、应力集中[6]和平均应力[7]等均对疲劳寿命产生重要影响．应力梯度对缺口件疲劳寿命有明显的影响，当缺口处有较大的应力梯度时，如简单地使用危险点处的应力应变进行寿命预估时其结果会相对保守[8]．这是由于缺口处应力场分布不均,当结构表面高应力区材料达到屈服强度且产生塑性变形时，邻近区域的材料仍处于较低的应力水平，对高应力区材料起到一定的支撑作用，从而减缓了裂纹的萌生、扩展，以及最终导致疲劳失效的进程[9]．文献[6]根据临界距离的思想，考虑了多轴加载时应力梯度效应对材料疲劳损伤演化过程的影响，提高了多轴缺口件疲劳寿命的预估精度．文献[10]在研究中周疲劳过程中，将临界距离视为与疲劳寿命相关的函数．文献[11]在文献[10]的基础上引入应力集中系数对该函数予以修正．上述模型在对疲劳寿命进行评估过程中不能事先预知临界距离的具体数值．文献[12]认为临界距离是材料常数，其数值与缺口形状、疲劳寿命及所施加的外部载荷均无关．文献[13]认为临界距离可以表述为与载荷比和峰值应力相关的函数，与疲劳寿命无关．目前，临界距离理论在进行缺口疲劳分析中有较强的实用性，但临界距离的确定仍是一个难题[9]．能量法认为外载荷每循环一次能量就在构件里累积一次，从而对构件造成不可逆转的损伤，当损伤累积到一定程度时，材料发生疲劳破坏[14]．然而当构件服役时，并非所有材料平面的损伤都是相同的，其损伤程度与该平面上的应力应变分量有关，当某一平面内所吸收的应变能最大时，则裂纹沿此方向萌生的概率最大[15]．以能量为疲劳损伤的参量，以临界面作为裂纹萌生的方向，采用较少的待定参量，对疲劳损伤参量的定义更加明确，克服传统能量法作为标量难以描述裂纹扩展方向的缺点，具有明确的物理意义．本研究以多轴比例加载条件下的缺口件为研究对象，基于临界距离的思想，考虑应力梯度对构件疲劳损伤造成的影响，给出一种确定疲劳损伤影响区的理论计算方法，结合临界面理论，以裂纹萌生面上疲劳损伤影响区的场强为损伤参量进行寿命预估，并与试验结果和局部应力应变法进行对比，验证本文方法的可行性．1 多轴缺口件临界面的确定首先利用弹塑性有限元软件ANSYS对构件进行力学分析，根据材料硬化规律，采用随动硬化模型描述材料的包辛格效应，通过载荷步设置实现循环加载，获得单个循环内缺口根部应力应变历程；接着确定构件危险点的位置，提取危险点在基本坐标系xyz下的应力应变分量，通过坐标变换矩阵得到过危险点任意平面的应力应变分量；最后利用能量法给出任意平面的塑形应变能函数f(θ,φ)，利用该应变能函数分别对θ和φ求导得到其驻点，由θ和φ值最终确定临界面的位置．1.1　临界面坐标变换为了得到危险点在任意平面上的应力应变分量，可利用坐标变换的方法：将x-y平面绕z轴旋转角度θ得到新的坐标平面x'-y'，将z轴绕x'旋转角度φ得到新的坐标轴z'，结果如图1所示．.F001图1坐标变换示意图基本坐标系xyz下的应力张量σ和应变张量ε分别为：σ=σ11σ12σ13σ21σ22σ23σ31σ32σ33;ε=ε11ε12ε13ε21ε22ε23ε31ε32ε33．坐标变换矩阵为M=cos θ-sin θcos φsin θsin φsin θcos θcos φ-cos θsin φ0sin φsin φ．通过坐标旋转矩阵变换可以得到新坐标系x'y'z'下过危险点任意平面的应力张量σ'和应变张量ε'分别为σ'=MTσM；ε'=MTεM．能量法一般以能量作为损伤参量，但能量为标量，不能解释裂纹沿临界面萌生和扩展这一现象．而坐标变换矩阵可以解释能量法中疲劳裂纹萌生和扩展的物理意义，为能量法确定临界面位置提供理论依据．1.2　临界面的确定在多轴加载条件下，缺口根部典型单元应力状态如图2所示．.F002图2任意平面的应力分量任意平面MNOP上的应力分量为σ11'，τ12'，τ13'，其应变能密度关于θ和φ函数为f(θ,φ)=(σ11'ε11'+σ12'ε12'+σ13'ε13')/2，式中：σ1i'=(σ11M11T+σ21M12T+σ31M13T)M1i+(σ12M11T+σ22M12T+σ32M13T)M2i+(σ13M11T+σ23M12T+σ33M13T)M3i;ε1j'=(ε11M11T+ε21M12T+ε31M13T)M1j+(ε12M11T+ε22M12T+ε32M13T)M2j+(ε13M11T+ε23M12T+ε33M13T)M3j,其中i，j分别为1，2，3．为简化过程，令(•)d=[(•)11-(•)22] /2，(•)m=[(•)11+(•)22] /2，•表示应力或者应变．此时应变能密度函数为f(θ,φ)=(1/2)[σmεm+σdεd+σ12ε12-σ33ε33-(σmεd+σdεm)cos(2θ)-(σmε12+σ12εm)sin(2θ)]sin2 φ+σ33ε33 . (1)鉴于f(θ,φ)为应变能密度函数，其极值点所在平面的应变能密度最大，疲劳裂纹也易在此处萌生．对式(1)求驻点，可得到f(θ,φ)极值点时的θ，φ值．具体求解过程如下．对φ求偏导可得∂f(θ,φ)/∂φ=[σmεm+σdεd+σ12ε12-σ33ε33-(σmεd+σdεm)cos(2θ)-(σmε12+σ12εm)sin(2θ)]sin φcos φ  .令∂f(θ,φ)/∂φ，得φ=90°或φ=0°．对θ求偏导可得     ∂f(θ,ϕ)/∂θ=0=[(σmεd+σdεm)sin(2θ)-(σmε12+σ12εm)cos(2θ)]sin2 φ .令∂f(θ,ϕ)/∂θ=0可得：a. φ=mπ,m=0,1,2，θ为任意值，此种情况与实际不符，舍弃；b. φ≠ mπ时，tan(2θ)=(σmε12+σ12εm)/(σmεd+σdεm)，与实际情况相符，因此临界面位置为θ=12arctanσmε12+σ12εmσmεd+σdεm,   φ=90°．(2)借助ANSYS提取危险点处的应力应变分量，结合坐标变换矩阵和式(2)即可计算临界面的位置．2 考虑应力梯度的多轴疲劳损伤参量无论构件承受静载荷还是循环载荷，材料的结构强度由缺口强度决定，由于缺口效应的存在，缺口危险点周围的应力应变场存在应力应变梯度，应力应变梯度对材料强度和疲劳寿命均有影响[16]．姚卫星等[17]依照材料损伤机理和疲劳损伤的细微观与宏观分析，提出了场强法来反映缺口件受载的损伤情况，并假设光滑试件的应力应变场强的历程和缺口根部的应力应变场强历程相同，则二者的疲劳寿命相同，具体为σF1=1V∫Ωf(σij)ϕxdν，(3)式中：σF1为应力场强度；Ω为疲劳破坏区；V为Ω的体积；f(σij)为破坏函数；ϕx为权函数．以临界面上等效应力为基础考虑应力梯度效应的影响，因此式(3)可以改写为σF1=1S∫Df(σeq)ϕxdx，(4)式中：S为临界面上疲劳损伤区D的面积；f(σeq)为临界面上等效应力分布函数．对于二维结构，权函数可以表示为ϕx=1-x|χx|(1+sin θ)  ,式中：x为积分点到危险点距离；θ为临界面位置与试件径向的夹角，这里定义临界面为等效应力梯度积分方向；χx=f-1(σeq)df(σeq)/dx为积分点处相对应力梯度．对于缺口件而言，在过危险点的剖面上危险点处应力最大，之后随着x的增加而减小，从而使构件实际寿命与利用危险点峰值应力作为损伤参量得到的结果不一致[18]．为克服这一不足，将疲劳损伤区的场强作为疲劳损伤参量，该参量考虑缺口附近应力梯度对疲劳寿命的影响．为进一步简化计算过程，以临界面上某一特定方向为等效应力场强积分方向，可以将式(4)的面积分转化为线积分形式σes=1R∫0Rf(σeq)ϕxdx，(5)式中：σes为疲劳损伤区的等效应力场强；R为临界面上疲劳损伤区内给定路径的临界距离值．为确定R的大小，对相对应力梯度函数求驻点，即dχx/dx=0．(6)利用式(6)得到相对应力梯度取极值时距危险点处的距离xeff，并定义此时的xeff=R，再借助式(5)计算多轴比例加载下缺口件疲劳损伤影响区的场强，结合Manson-coffin方程可以预估疲劳寿命．3 多轴疲劳试验试验材料选用Q345钢，其弹性模量E=1.92×1011 Pa；剪切模量G=7.38×1010 Pa；屈服应力σy=476 MPa；抗拉强度σu=625 MPa；泊松比ν=0.3；伸长率δ=22%．该材料单轴疲劳性能参数为：循环强化系数K=2.85×109 Pa；疲劳强度系数αf'=1.44×109 Pa；循环硬化指数n=0.378；疲劳强度指数b=-0.159；疲劳延性系数αf=0.249；疲劳延性指数c=-0.494．试验件的基本形状参考GB3075—82，为考虑缺口几何参数对多轴疲劳寿命的影响，同时结合试验机夹具的要求，缺口试件选用U型缺口(图3)，每个工况试验三次取其平均值，试验结果如表1所示，表中：r为缺口半径；t为缺口深度；εa为轴向应变；εb为切向应变．疲劳试验在Instron 8850拉扭疲劳试验机上进行，利用拉扭引伸计Epsilon 3550控制轴向和切向应变，试验过程中轴向应变和切向应变的相位差为0°(比例加载)，拉压和扭转加载频率均为1 Hz．试件断裂认为疲劳破坏产生，此时的循环周次记录为J．.F003图3缺口试验件(mm).T001表1多轴比例疲劳试验结果r/mmt/mm工况εa/%εb/%J1.603.6L10.3950.47912 0671.603.2L20.3950.47912 8651.602.4L30.3950.47913 6801.601.6L40.3950.47914 9361.284.0L50.3950.4798 9631.922.08L60.3950.47915 6301.283.6L70.3470.40812 9341.283.6L80.3470.34713 3911.283.6L90.3470.29515 5631.923.6L100.3470.40812 1521.923.6L110.3470.24513 8611.923.6L120.3470.20012 9622.243.6L130.3470.40812 0302.243.6L140.2950.40814 1652.243.6L150.2510.25113 3944 多轴缺口件应力应变分析典型算例以L1工况作为典型算例对本文方法进行详细描述，其他工况仅给出具体结果．利用ANSYS对试件进行应力应变计算时，采用多线性随动强化模型MKIN．划分网格时采用自由网格划分方法，缺口处自动加密．设置边界条件和施加载荷时，固定试件一端，在另一端分别施加周向载荷和轴向载荷，其中周向载荷转化为扭转角，轴向载荷转化为位移均匀施加在右端面最外层圆周的所有节点上．L1工况下等效应力云图如图4所示．图中等效应力最大值625 MPa出现在加载端，这是载荷施加在节点上导致的，依据圣维南原理，不影响缺口附近的结果．缺口根部最大等效应力σp=564    MPa并定义为危险点，其中临界面的位向如图4中AC所示，A为危险点位置．.F004图4L1缺口试样的等效应力云图5 结果与讨论为分析L1工况下临界面上等效应力分布，在ANSYS中将工作平面旋转θ=34°，使工作平面与临界面重合，沿工作平面将试件进行剖分，给出临界面上等效应力分布云图(图5)．.F005图5L1工况下临界面上应力分布以危险点A为起点，沿图5中AB路径提取应力值，路径长度为5 mm，借助Matlab软件对提取的等效应力值进行函数拟合，拟合时只须关注距缺口根部一定范围内拟合函数与有限元提取结果的符合程度，即只要在一定范围内拟合函数能够表征应力的分布即可．这里采用6次多项式进行拟合：f(σeq)=B1x6+B2x5+B3x4+B4x3+B5x2+B6x+B7. (7)L1工况下从软件中提取B1~B7分别为-0.074 2，1.408 0，-10.458 5，35.204 3，-28.461 2， -139.080 9，565.427 7．有限元计算结果和拟合函数如图6所示，可以看出：在距缺口根部5 mm内，拟合曲线可以很好地表征等效应力的分布状态，临界面上等效应力梯度为df(σeq)/dx=6B1x5+5B2x4+4B3x3+3B4x2+2B5x+B6.相对等效应力梯度为χx=f-1(σeq)df(σeq)/dx．(8)L1工况下使式(8)取极小值的xeff=0.71，即R为0.71 mm(图6)．将等效应力分布函数(式(7))、θ和R值代入式(5)中即可以计算疲劳损伤区的等效应力场强σ es，结果如表3所示．多轴比例加载下应力应变关系可以用Osgood-Ramberg方程表述，具体为Δε/2=Δσ/(2E)+[Δσ/(2K)]1/n，式中：Δε/2为应变幅值；Δσ/2为应力幅值；n为循环应变硬化指数；K为循环强度系数．.F006图6临界面上相对应力梯度借助Osgood-Ramberg方程可以给出临界面上等效应力场强-等效应变场强的关系式为εes/2=σes/(2E)+[σes/(2K)]1/n．(9)由式(9)可得到相对应的等效应变场强εes，结合Manson-coffin方程式εeq=(σf'/E)(2Nf)b+εf'(2Nf)c可以得到考虑应力梯度效应的多轴比例加载下缺口试件的疲劳寿命Spt．由缺口根部峰值应力σp，结合Osgood-Ramberg方程和Manson-coffin方程可以得到利用局部应力应变法计算的疲劳寿命Spl，结果如表3所示．.T002表3缺口件疲劳寿命工况R/mmθ/(°)σes/MPaσp/MPaSptSplL10.7134521.8564.213 1438 385L20.7231518.4558.313 6508 905L30.6932514.5551.914 2609 514L40.8037509.5547.615 0919 951L50.6832551.4604.89 5645 642L60.7535502.4548.616 3739 847L70.8230515.9559.914 0388 760L80.7638510.0553.615 0069 347L90.7239499.8542.516 87510 501L100.7038520.3564.713 3648 342L110.6931513.8558.814 3738 859L120.8137505.6549.915 7809 714L130.6937516.5559.513 9438 796L140.8232508.2551.615 3179 544L150.7835507.9552.415 3699 465从表3可以看出：在多轴比例加载情况下，考虑应力梯度后的场强明显低于危险点处的峰值应力，分别以疲劳损伤影响区的场强和危险点处的峰值应力作为疲劳损伤参量对疲劳试件进行寿命预估，并与试验结果进行对比，结果如图7所示．图7表明：本文方法与局部应力应变法预估寿命的误差均在2倍因子内，但本文方法预估精度要高于局部应力应变法，这是因为在缺口附近存在明显的应力梯度现象，内部低应力区的晶粒对外层高应力区的晶粒有一定的支撑作用，进而影响危险点处的疲劳损伤演化过程．而局部应力应变法是以危险点处的峰值应力作为损伤参量进行寿命预估，没有考虑缺口附近区域应力梯度效应的影响，导致预估精度相对偏低，且结果偏保守．.F007图7预估结果与实验结果对比