金属梁是工程中应用广泛的重要结构件，在复杂环境中长期经受周而复始的弯曲变形．在高频次(高周)与低频次(低周)反复弯曲过程中，梁作为承受弯曲载荷构件，其抗弯强度将逐渐弱化直至最终丧失而发生失效破坏．然而，在每次弯曲循环中，该强度弱化效应极其微小，通常难以察觉，造成事故突发．由于在反复变形过程中，金属材料强度的渐变弱化是导致工程结构失效破坏的主因[1-3]，因此探讨金属梁的弯曲疲劳破坏特征及其规律具有重要意义．金属梁疲劳失效问题涉及金属材料弹塑性变形特征、截面几何特征、应力与应变的非均匀分布及载荷的循环变化特征等多种因素，根本性因素是由于金属材料不可避免的耗散累积效应导致其强度在循环过程中渐变弱化．事实上，该因素导致下述复杂性：首先，梁截面上非均匀应力分布特征并非随着循环进程保持不变，而是在每次循环后均将发生变化；其次，由于该变动，弯矩-曲率关系随循环进程也将发生变化，从而导致梁弯曲强度在循环进程中逐渐下降直至最终丧失．如上所述，循环弯曲金属梁的疲劳破坏模拟本质上在于如何准确表征其弯曲强度在循环弯曲过程中的渐变弱化直至最终丧失．通常方法往往基于金属材料理想弹塑性特征的粗略近似，完全略去金属的强化和软化效应．由此只能导出一些初步结果，无法模拟典型循环弯曲过程(截面应力远低于屈服强度)中梁弯曲强度的渐变弱化效应．在金属的疲劳失效表征方面，遵循现有途径不仅须分别处理高中低周各种情形，而且只能给出适于某些特定情形的结果．现有途径主要基于断裂力学方法[4-11]和损伤力学方法[12-15]，并引入特设失效判据．基于损伤变量可表征强度的渐变弱化，但须寻求损伤变量演化方程的恰当形式及决定大量关联于损伤的未知参数．关于疲劳失效的模拟方法及其特设判据等结果参见文献[6，16-19]．基于现有途径难以表征梁截面上非均匀应力分布在循环弯曲过程中的复杂变化及相关疲劳失效问题，须借助有限元手段，然而经由基于有限元步骤的反复计算将耗费大量的时间．基于全程耗散新弹塑性模型[20-29]，本研究给出金属梁高低周弯曲疲劳破坏的直接模拟．在既不涉及通常的损伤内变量又不涉及裂纹演化及相关特设失效判据的基础上，直接表征金属梁在循环弯曲过程中其弯曲强度全程渐变弱化直至丧失的复杂行为．提出高效算法导出金属梁在循环弯曲过程中截面上非均匀应力-应变和弯矩-曲率的显式关系，并分析金属梁的非线性循环弯曲响应及强度弱化效应，由此直接导出破坏循环数随曲率幅变化的疲劳特征曲线并与实验数据进行对比．1 全程耗散新弹塑性模型引入欧拉率型弹塑性本构表述，以变形率D的加法分解为出发点[30-31]D=De+Dp，(1)式中：De为弹性变形率；Dp为塑性变形率．1.1　自洽弹性率方程分解式(1)中的弹性变形率由下述弹性率方程决定[32]De=12Gτ˙log-νE(trτ˙log)I，(2)式中：G，ν和E分别为弹性剪切模量、泊松比和弹性模量，且有E=2G(1+ν)；I为二阶单位张量；τ˙log为引入克希霍夫(Kirchhoff)应力张量τ的对数共旋率[32]，τ˙log=τ˙+τ⋅Ωlog-Ωlog⋅τ，(3)Ω为共旋率．已证明[32]，弹性率方程(2)在下述意义上自洽：在纯弹性情形下该率方程精确，可给出基于对数应变的有限超弹性应力-应变关系．1.2　塑性率方程分解式(1)中的塑性变形率，由下述正交流动法则[33]指定Dp=ρζ∂f∂τ，(4)式中f为von Mises函数，f=g-r, g=tr(τ˜2/2), r=q2/3，(5)其中，τ˜为偏量应力，q=q(k)为随塑性功k演化的屈服强度，k˙=τ˜:Dp．(6)式(4)中的ζ包括加载项f^和塑性强度项h^，塑性因子ρ决定塑性流，并决定耗散效应．1.3　连续塑性因子与全程耗散按照经典弹塑性表述，塑性因子ρ只能取离散值0和1，即在加载情形取值1，而在卸载情形取值0．加载和卸载情形基于屈服条件与应力增量取向予以判断(参见文献[34])．只有当应力点保持在屈服面上时才会产生塑性流，而当应力点位于屈服面内时则不出现塑性流．研究指出[20-23,27]：上述离散塑性因子导致从弹性态转变到塑性态时发生模量强间断．或许更为本质的是，任何基于屈服的弹塑性模型排除了循环加载下的疲劳失效效应[20-22,25-27]，而疲劳失效效应正是金属材料的典型非弹性特征．实际上，对金属的实际变形行为，即使远低于初始屈服强度的应力水平也将诱导塑性变形(参见文献[35])，对应于0~1的塑性因子．引入光滑取值的塑性因子ρ=emg/r-1em-1，(7)式中m为失效循环次数．式(7)给出的新塑性因子在应力点接近于经典屈服面时近似于1，而当应力点远离经典屈服面(即g/r≪1)时趋于零，当应力点位于经典屈服面上时恰为1，而应力为零时为零．基于式(7)给出的光滑塑性因子，可建议新塑性率方程如下[20-29]：Dp=ρf^+|f^|2h^∂f∂τ;       f^=2Gτ˜:D;h^=4g(3G+qq')/3. (8)上述方程不再如经典表述那样将基于屈服的加卸载条件引入为外加强制性条件，而是将这些条件自动包含为固有响应特征[20]．依照新塑性率方程(7)和(8)，当应力点接近于经典屈服面时塑性变形将较为显著，而当应力点远离经典屈服面时将非常小．方程(7)和(8)同时将疲劳失效效应自动纳入为下述渐近条件下的固有响应特征[20,22]limk→∞q=0，(9)该条件意味着塑性强度将随塑性功(耗散)累积而渐近弱化直至丧失．联合式(1)、(2)、(5)、(6)和(8)，得到新弹塑性模型如下：D=12Gτ˙log-νE(trτ˙)I+ρf^+|f^|2h^∂f∂τ;                    k˙=ρf^+|f^|2h^τ˜:τ˜. (10)不同于基于屈服的经典弹塑性模型，新模型在应力幅低于屈服强度的每一循环加载过程中将给出耗散单调累积效应，于是在循环过程中耗散单调增长而将自动达到下述临界失效条件[22-29]3G(1-ρ)+qq'=0，(11)因而对每一循环过程，全程耗散新弹塑性模型将自动预言疲劳失效效应．为此，只须适当选取满足渐近条件式(9)的应力极限q，然后由式(11)判断何时发生失效破坏．考虑下述简单形式q=12q01-tanh βkkc-1，(12)式中q0，kc和β分别为应力极限的最大值、塑性功(耗散)的临界值及渐近软化效应．无量纲正参数β0表征软化过程中的延性效应，高延性(平缓应力-应变曲线)对应较小β值．kc为塑性功k的临界值，当k越过kc后，q将快速趋于零．2 单调单轴响应与参数识别式(12)中的参数由单调单轴响应数据决定，这将是后续分析的出发点．设h和τ是金属杆的轴向对数应变和轴向克希霍夫应力，对单调单轴变形响应有τ=τ(h)，k=k(h)，于是上述响应可由方程(10)的一维简化形式即下述微分方程组决定：dhdτ=1E+ρ3G(1-ρ)+qq';     dkdτ=ρτ3G(1-ρ)+qq';        τ0=0,  k0=0. (13)给定包括软化效应的一组单调单轴实验数据(hs，τs)，s=1，2，…，N，可找到参数m，q0，kc和β的适当值，使得由方程组(13)决定的应力-应变响应拟合该组数据．不同于通常只涉及强化效应的参数识别，为达到预期目标，须包含软化效应的实验数据．金属单轴试样在测试中将出现颈缩现象，导致非均匀应力应变分布，因而难以获得软化数据．文献[36]提出了处理该困难问题的解决方案，依据该方案可获得直至失效的完整单调单轴实验数据．获得适当参数值后，即可基于新模型模拟循环弯曲响应．3 循环弯曲的表征考虑具有对称截面的均匀金属梁的纯弯曲变形．设2H为梁高，高度坐标为z处梁的半宽度y= y(z)，-H≤z≤H．在每次循环中弯曲曲率α从某个值A*增大到给定幅值A，此处A*包括两种情形，其一为对称循环弯曲，即A*=-A；其二为自然循环弯曲，即每次循环中达到给定幅值A后将所加弯矩卸除为零，于是A*为每次循环后产生的塑性曲率．设曲率值为A0时梁截面上最大法应力刚好达到经典意义上的初始屈服应力．对通常的循环弯曲过程，曲率幅值A远低于临界值A0，于是基于屈服的经典弹塑性模型，在循环弯曲过程中梁截面上法应力分布τ(z)应当总是线弹性的且对所有循环均相同，从而排除了任何疲劳失效效应．若不引入损伤内变量和裂纹演化以及特设失效判据等附加假设，则将无法模拟疲劳失效效应．另一方面，依照全程耗散新弹塑性模型，即使是极低的应力水平也可产生微小的塑性应变，同时塑性功(即耗散)k将随循环数增加而单调累积增大．由此可见，不同于基于屈服的经典弹塑性模型，新模型预言：梁截面上法应力分布τ(z)在每个循环中是非弹性的，且随着循环进程而演变，导致梁截面抗弯刚度随着循环进程而渐变弱化直至疲劳失效．截面高度坐标为z处每点的法应变为ε(z)=αz    (A*≤α≤A,-H≤z≤H)，(14)式中α为弯曲曲率．对小曲率弯曲，上述应变给出对数应变．高度坐标为z处每点的法应力与塑性功的控制方程为：  d[τ(z)]d[ε(z)]=1E+ρ(z)3G[1-ρ(z)]+ξ(z)-1;d[k(z)]d[ε(z)]=d[τ(z)]d[ε(z)]ρ(z)τ3G[1-ρ(z)]+ξ(z);                      A*z≤ε(z)≤Az;                     ξ(z)=q(z)q'(z);   q(z)=12q01-tanhβk(z)kc-1;              ρ(z)=exp[mτ2(z)/q2(z)]-1em-1. (15)4 循环弯曲响应与疲劳特征曲线基于方程组(14)和(15)可直接模拟梁循环弯曲的全程响应．对每一给定高度坐标z，由于累积效应，塑性功在每个循环起始时的初值由其在上一循环的终值给出．于是，从方程组(15)导出的截面应力分布τ(z)将随循环进程而演变，由此即可模拟梁弯曲刚度随循环进程的渐变弱化．在每次循环中对应于曲率幅值A的最大弯矩将随循环进程逐渐下降，这对应于梁弯曲强度在循环进程中的渐变弱化．对给定曲率幅值A，设N是达到失效时的循环数，此即对应于曲率幅值A的疲劳寿命．导出疲劳寿命N随曲率幅值A的变化曲线有重要意义．4.1　应力分布的高效算法设kn(z)和kn+1(z)是第n次循环的初值和终值．注意到每次循环中塑性功的增量非常小，应用新近提出的直接递推算法[30]，可得到下述高效计算方案kn+1(z)=kn(z)+γΔkn+1(z);k0(z)=0  (n=0,  1,     …,  N), (16)式中γ为循环因子，对自然循环和对称循环分别取值1和2．τn+1=τn+1 (z)为第n+1次循环中对应于曲率幅A在高度z处的法应力，则有：           Δkn+1(z)=qn2(z)(3G-bn)6Gmbn∙    ln3G[1-ρn(z)]+ξn(z)3G+ξn(z)-τn+12(z)2bn;              bn=3Gem+(em-1)ξn(z);          ρn(z)=exp[mτn+12/qn2(z)]-1em-1;     qn(z)=12q01-tanhβkn(z)kc-1;                    ξn(z)=qn(z)∂q(z)∂kk=kn(z). (17)基于上述表达式，考虑方程(15)的三阶摄动解，可导出截面上法应力分布的显式表达式如下：           τn(z)=Eε(z)-Bn(z)ε3(z); Bn(z)=mE43qn2(z)[3G+ξn(z)](em-1). (18)4.2　弯曲强度的渐变弱化对第n次循环，由式(14)和(18)可导出应力-曲率关系如下τn(z)=Ezα-Bn(z)(zα)3．(19)梁弯曲强度的演化由弯矩(M)-曲率关系表征，对第n次循环，应用式(19)可导出该关系如下Mn(α)=4∫0H[Ezα-Bn(z)(zα)3]y(z)zdz．(20)式(20)表明：由于塑性功的累积效应，对应于曲率幅α=A的最大弯矩随循环数n变化且实际上随循环数n增大而下降．这意味着在循环进程中梁弯曲强度将随耗散累积而逐渐弱化．然而，对小曲率幅值A，在两个相继循环中所述变化极其微小．无论该变化多么微小，充分大循环数的耗散累积将导致疲劳失效．在给定曲率幅的循环弯曲情形，当最大弯矩历经循环数N后下降到某一给定百分数(如50%)时，可认为失效即将发生．设金属梁的截面形状和构成材料给定，则对每个给定曲率幅A的循环弯曲过程，可由式(20)计算该金属梁的弯矩-曲率响应，特别地，可计算直到失效时的循环数N．对一系列曲率幅的值，如此将决定对应的一系列失效循环数，从而获得梁弯曲循环的疲劳特征曲线．5 数值算例5.1　单调单轴实验数据的拟合与参数识别如前所述，相关材料参数由拟合单调单轴实验数据决定．下面分别考虑由Dowling等[37]给出的关于RQC-100钢单轴数据及由Zhou等[38]给出的贝氏体钢单轴数据，参数值由表1给出，单轴实验数据如图1所示，图中：ε为应变；τ为应力．由图1可知：拟合结果符合显著强化前的初始数据及强化后的软化数据．式(12)给出的应力极限表达式中不包括显著强化效应，出现与强化阶段数据的偏离在引入强化项后可得到进一步的拟合结果[27-29]．10.13245/j.hust.239334.T001表1拟合单轴数据的参数值材料E/GPaνq0 /MPakc/MPamβRQC-10钢2000.38575 4500.70.8贝氏体钢1600.31 6202950.16.510.13245/j.hust.239334.F001图1RQC-100钢和贝氏体钢的单轴实验数据5.2　RQC-100钢梁的对称循环弯曲考虑RQC-100钢梁(单轴实验数据参见文献[37])，矩形截面尺寸为12.7 mm×6.3 mm．应用上节高效算法计算在对称循环弯曲下的模拟结果，图2为弯曲疲劳特征曲线．在A=1.58，0.96，0.69 m-1下最大弯矩随循环数的变化曲线与实验数据对比如图3所示，图中M为最大弯矩．模拟结果均与Dowling等[37]给出的实验数据一致．其中，对较大曲率幅值即A=1.58 m-1，当循环数较大时模拟结果与实验数据出现一定程度偏离．这主要源于决定第n次循环中截面应力分布的式(19)略去了高阶项．因而当曲率幅值较大时，误差随循环数增长而增长．预计包含高阶项的新表达式可减少此误差．10.13245/j.hust.239334.F002图2RQC-100钢梁的弯曲疲劳特征曲线10.13245/j.hust.239334.F003图3最大弯矩随循环数的变化曲线与实验数据的对比5.3　贝氏体钢梁的自然循环弯曲考虑贝氏体钢方形梁(单轴实验数据参见[38])，方形截面边长为500 mm．处理自然弯曲循环，即每次循环中达到给定曲率幅值A后将所加弯矩卸除为零．应用上节算法计算在自然循环弯曲下模型的结果．图4给出A=0.000 3，0.000 6，0.000 9 m-1下直至疲劳失效的全程弯矩-曲率响应．10.13245/j.hust.239334.F004图4贝氏体钢梁在自然循环弯曲下全程弯矩-曲率响应如图4所示，每次循环的最大弯矩随循环进程逐渐下降．当达到某个循环数，最大弯矩下降到其初值的某个指定百分数(取为50%)时，疲劳失效实际上即将发生．6 结语基于无屈服新弹塑性模型导出了循环弯曲全程响应的控制方程组，给出了直接积分该方程组的高效数值算法，导出金属梁在循环弯曲过程中截面上非均匀应力-应变和弯矩-曲率的显式关系，并分析金属梁的非线性循环弯曲响应及强度弱化，由此直接导出破坏循环数随曲率幅变化的弯曲疲劳特征曲线及弯矩-曲率响应结果，无须借助较为耗时的有限元手段．数值结果表明，模型预言与梁循环弯曲的实际疲劳特征一致．后续研究须考虑金属梁循环弯曲的其他方面，包括弯矩循环下疲劳失效特征及非矩形截面等．此外，须进一步考虑钢筋混凝土梁，特别是拉压非对称变形行为及中性轴在循环进程中的移动等．