一种基于损伤力学非概率区间分析模型的金属疲劳裂纹全寿命预估方法

摘要

本发明公开了一种基于损伤力学非概率区间分析模型的金属疲劳裂纹全寿命预估方法。该方法首先在损伤力学模型中选择一种损伤演化方程建立结构包含损伤信息的有限元分析列式，查询疲劳试验手册拟合得到损伤演化方程中的参数；然后结合损伤力学有限元法与区间有限元法，将初始损伤度与损伤参数看作区间不确定变量以表征疲劳寿命分散性；进一步建立结构的有限元分析模型，给定初始临界单元损伤度增量，不断迭代增加所有单元的损伤度，通过各个单元损伤度的大小判定单元破坏并重新赋予强度与刚度属性；最后当裂纹扩展达到临界裂纹长度后判定结构破坏，由损伤演化方程变式与区间顶点传播分析方法计算得出疲劳寿命范围。本发明更加精细化的预估了疲劳裂纹寿命。

说明书

技术领域

本发明涉及疲劳断裂和损伤力学领域，特别涉及考虑不确定性作用下裂纹扩展尺度对结构安全性能的定量表征以及基于损伤力学有限元与非概率区间有限元结合下的疲劳裂纹全寿命精细化分析。

背景技术

随着科技的发展，结构的设计思想已经发生经历了静强度、动强度、疲劳强度和断裂强度的演变。自损伤力学产生以来，用损伤力学方法研究疲劳裂纹问题成为一种新的趋势。损伤力学是一门较系统的研究微缺陷以及这些缺陷对应力和应变状态影响的科学，是固体力学研究领域的一门新兴学科，其理论基础是固体力学和不可逆过程的热力学。这两个理论能够成功地用来解释材料性能，而无需详细讨论材料物理微结构的复杂性。预测构件疲劳损伤寿命的损伤力学计算方法从本质上分析疲劳裂纹的形成和扩展，发展和完善了疲劳断裂理论，在实际工程中应用很方便。

然而，构件或者材料的疲劳裂纹扩展行为研究涉及了力学、材料、机械设计与加工工艺等诸多学科，影响疲劳裂纹扩展的因素也非常之多，包括裂纹的几何形态、初始裂纹长度、材料特性、裂纹扩展规律、扩展方向、构件的几何尺寸和载荷历程等。由于初始参数的分散性，应用数值方法预估疲劳裂纹扩展寿命的结果必然会有误差。建立一个包含上述各影响因素的疲劳裂纹分析模型，并且准确预测疲劳裂纹全寿命是一件困难的事情。工程结构的服役环境相对复杂，制造加工工艺及材料非均质性所造成的初始缺陷和损伤不可避免，并在未来长期服役过程中于结构内部不断发展、蔓延、传播，严重影响着结构的力学行为及使用安全。综合上述情况，为分析数值方法预测的精度，就很有必要精细化研究预测裂纹的全寿命。

当前，随机建模及数值方法在不确定性分析领域发挥了重要作用，但用随机理论求解问题时，事先需要知道大量的试验信息确定模型输入参数的概率分布规律。在实际工程中，获取充足的试验数据往往代价昂贵。如此一来，信息的缺乏使得概率模型不能真实的描述客观实际，这在一定程度上限制了随机模型的应用。

因此，使用非概率区间分析方法来表征参数的不确定信息，基于初始参数的不确定性，研究不确定性传播问题导致预测裂纹全寿命的影响程度，精细化裂纹全寿命预估方法具有显著的现实意义。

发明内容

本发明要解决的技术问题是：克服现有技术的不足，提供一种基于损伤力学非概率区间分析模型的金属疲劳裂纹全寿命预估方法。本发明充分考虑实际工程问题中普遍存在的不确定性因素，以提出的非概率区间顶点方法分析不确定性传播问题，所得到的裂纹全寿命范围符合真实情况，工程适用性更强。

本发明采用的技术方案为：一种基于损伤力学非概率区间分析模型的金属疲劳裂纹全寿命预估方法，实现步骤如下：

第一步：根据金属材料类别与构件所承受载荷形式选择一种损伤演化模型，单轴加载条件下损伤演化方程可表示为：

$\frac{dD}{dN}=\frac{\alpha }{p+1}{\left(\frac{1}{2E}\right)}^{p+1}{\left(\frac{1}{1-D}\right)}^{2p+2}[{(1-R)}^{q(p+1)}{\sigma }_{Me}^{2p+2}-{\sigma }_{th0}^{2p+2}{(1-D)}^{(0.5+\beta )(2p+2)}]$

其中，D代表在0与1之间变化的单元标量损伤度，N代表应力循环次数即单元寿命，β、α和p代表材料的损伤力学参数，可由材料的疲劳性能曲线确定，E代表材料弹性模量，R代表循环载荷的应力比，σ_Me代表材料受到最大载荷时单元对应的等效应力，σ_th0为无初始损伤情况对应的应力门槛值；

第二步：将损伤演化方程耦合传统有限元方法，得到给定损伤场时结构应力分析的损伤力学有限元分析列式：

$\left(\underset{e}{\Sigma }{\left[A_e\right]}^T\right[K_e\left]\right[A_e]-\underset{e}{\Sigma }{D}_e{\left[A_e\right]}^T[K_e\left]\right[A_e\left]\right)\left\{\delta \right\}=\left\{f\right\}$

其中，[A_e]为位移协调矩阵，[K_e]为单元刚度矩阵，D_e为单元损伤度，{δ}为总体位移列阵，{f}为载荷列向量；

第三步：获取标准疲劳试验件的中值疲劳寿命数据拟合损伤演化方程中的损伤参量，将损伤演化方程在0到1上积分，对应一个单元破坏则是标准的S-N曲线的寿命值，然后用最小二乘法拟合β、α和p；

第四步：利用区间向量x∈x^I＝(△D,β,p,α)合理表征贫信息、少数据条件下的结构不确定性，这里△D代表临界单元的损伤度增量，于是有：

x^U＝(△D^U,β^U,p^U,α^U)＝(△D^c+△D^r,β^c+β^r,p^c+p^r,α^c+α^r)

x^L＝(△D^L,β^L,p^L,α^L)＝(△D^c-△D^r,β^c-β^r,p^c-p^r,α^c-α^r)

其中，损伤力学参数β、α和p分别表示为区间变量，上标U代表参量的取值上界，上标L代表参量的取值下界，上标c代表中心值，上标r代表半径；

第五步：建立有限元模型，施加边界条件，先计算初始损伤场均为零，即无损伤情况下的应力场，利用二次开发编写程序提取各个单元等效应力存储于数组；

第六步：分析上一步提取得到的单元应力数组，由相对损伤度的最大值判断选择临界单元，表示如下：

$\underset{i\in [1,n]}{max}[{\left(\frac{dD}{dN}\right)}_i/(1-D_i)]$

其中，表示单元的绝对损伤演化率；

第七步：引入区间传播分析的顶点法，选择不确定参数的顶点上下限进行非概率不确定性传播分析，顶点法可表示为：

$\left(\begin{array}{c}{U}_1=(\underset{‾}{\Delta D},\underset{‾}{\beta },\underset{‾}{\alpha },\underset{‾}{p})\\ U_2=(\underset{‾}{\Delta D},\underset{‾}{\beta },\underset{‾}{\alpha },\bar{p})\\ U_3=(\underset{‾}{\Delta D},\underset{‾}{\beta },\bar{\alpha },\underset{‾}{p})\\ .\\ .\\ .\\ U_{16}=(\overline{\Delta D},\bar{\beta },\bar{\alpha },\bar{p})\end{array}\right)$

其中，U₁…U₁₆代表不确定参数组合取值形式，△D，β，α，p分别代表输入参数的下限，分别代表输入参数的上限；

第八步：给定临界单元损伤度增量△D，由损伤演化方程与临界损伤度增量计算所有单元的损伤度增量，将与前一步损伤场叠加得到的单元损伤场施加于有限元模型中进行带有损伤的应力分析，再提取单元等效应力存储于数组，并不断重复第六步直到判断临界单元的损伤度到1时即认为该单元破坏，并提取破坏单元长度与单元破坏寿命，长度即单元边长，任一单元损伤度增量与寿命增量可表示为：

$\Delta D\left(x\right)=\Delta D\left(e_i\right){\left[\frac{1-D\left(e_i\right)}{1-D\left(x\right)}\right]}^{2p+2}\left\{\frac{{\sigma }_{Me}{\left(x\right)}^{2p+2}-{\sigma }_{th0}^{2p+2}{[1-D\left(x\right)]}^{(0.5+\beta )(2p+2)}}{{\sigma }_{Me}{\left(e_i\right)}^{2p+2}-{\sigma }_{th0}^{2p+2}{[1-D\left(e_i\right)]}^{(0.5+\beta )(2p+2)}}\right\}$

$\Delta N=\Delta D/\left\{\frac{\alpha }{p+1}{\left(\frac{1}{2E}\right)}^{p+1}{\left(\frac{1}{1-D}\right)}^{2p+2}\right[{\sigma }_{Me}^{2p+2}-{\sigma }_{th0}^{2p+2}{(1-D)}^{(0.5+\beta )(2p+2)}\left]\right\}$

其中，△D(x)代表任一单元的损伤度增量，△D(e_i)代表临界单元的损伤度增量，D(x)代表任一单元的损伤度，D(e_i)代表临界单元的损伤度，σ_Me(e_i)与σ_Me(x)分别代表临界单元与任一单元的等效应力，σ_th0代表裂纹扩展门槛值，△N代表临界单元寿命增量，△D为临界单元的损伤度增量，E为弹性模量，α、β与p为损伤参量；

第九步：结合损伤力学有限元与区间顶点传播分析方法，将每一步迭代破坏单元的弹性模量将为极小值，并计算破坏单元的总长度a_i与材料的临界裂纹扩展长度a_c比较判断结构破坏，当a_i≤a_c即停止计算输出寿命的上下限与N。

其中，所述第一步中损伤演化方程的选择取决于结构几何、材料、载荷形式等输入参数的共同作用。

其中，所述第三步中通过查询疲劳试验手册获取标准疲劳试验件的中值疲劳寿命数据拟合损伤演化方程中的损伤参量，所述第三步中的损伤参数的拟合所需要的标准疲劳试验手册中值疲劳数据与构件真实载荷情况的应力比或者平均应力应该相同，对于手册中没有相同应力比应该通过应力幅与平均应力曲线进行等效转化，拟合数据时通过损伤演化方程与标准S-N曲线方程按最小二乘法拟合。

其中，所述第四步中区间不确定性参数向量x可以表示为：

x＝[x^L,x^U]＝[x^c-x^r,x^c+x^r]

＝x^c+x^r[-1,1]

＝x^c+x^r×e

其中，x^c×(ΔD^c,β^c,p^c,α^c)，x^r＝(ΔD^r,β^r,p^r,α^r)，e∈Ξ⁴，Ξ⁴定义为所有元素包含在[-1,1]内的4维向量集合，符号“×”定义为两个向量各对应元素相乘的算子，乘积仍为维数为4的向量。

其中，所述第五步中进行无损伤情况下的应力场分析，应该将所有单元的初始损伤度均设置为零，并将计算得到的单元Vonmises应力作为单元在外载荷下的最大等效应力。

其中，所述第五步中引入区间传播分析的顶点法，选择不确定参数的顶点上下限进行非概率不确定性传播分析，在引入顶点法进行传播分析时必须保证所研究的问题是单调的，针对裂纹萌生与扩展的全寿命分析随着迭代次数的增加损伤度和寿命均是单调递增的所以满足条件。

其中，所述第九步中临界裂纹扩展长度a_c为材料断裂韧性所决定的长度，是一个可查手册得到的确定值。

本发明与现有技术相比的优点在于：

(1)、本发明提供了裂纹扩展全寿命预估新思路，弥补和完善了传统分析方法的局限性。

(2)、本发明在进行预测过程中合理表征了初始损伤度与损伤演化方程材质参量的不确定性对金属材料疲劳裂纹全寿命的综合影响，并结合损伤力学有限元与区间有限元法得到了裂纹全寿命的区间范围，更加精细化的预估了疲劳裂纹寿命。

(3)、本发明使用的非概率区间顶点分析方法可以针对所有的损伤演化模型，相对于需要明确表达式传统的概率方法更为便捷。

(4)、本发明可对工程中关键结构件的疲劳裂纹全寿命进行预估，进一步可建立非概率可靠性模型，以此为基础建立优化方法指导结构设计可节约工程结构设计与试验成本。

附图说明

图1是本发明针对金属疲劳裂纹全寿命预估方法流程图；

图2是本发明针对损伤演化方程各参量拟合示意各参量拟合示意图；

图3是本发明标准试验模型尺寸示意图；

图4是本发明有限元网格与载荷约束示意图；

图5是本发明计算的裂纹扩展局部放大示意图；

图6是本发明针对LY12CZ板结构裂纹全寿命的上下范围与中心值示意图。

具体实施方式

下面结合附图以及具体实施例进一步说明本发明。

如图1所示，本发明提出了一种基于损伤力学非概率区间分析模型的金属疲劳裂纹全寿命预估方法，包括以下步骤：

(1)根据金属材料类别与构件所承受载荷形式选择一种损伤演化模型，针对经典的二维平板单轴加载条件下损伤演化方程可表示为：

$\frac{dD}{dN}=\frac{\alpha }{p+1}{\left(\frac{1}{2E}\right)}^{p+1}{\left(\frac{1}{1-D}\right)}^{2p+2}[{(1-R)}^{q(p+1)}{\sigma }_{Me}^{2p+2}-{\sigma }_{th0}^{2p+2}{(1-D)}^{(0.5+\beta )(2p+2)}]$

其中，D代表在0与1之间变化的单元标量损伤度，N代表应力循环次数即单元寿命，β、α和p代表材料的损伤力学参数，可由材料的疲劳性能曲线确定，E代表材料弹性模量，R代表循环载荷的应力比，σ_Me代表材料受到最大载荷时单元对应的等效应力，σ_th0为无初始损伤情况对应的应力门槛值。

对于不同的构件与载荷形式，损伤力学模型还有其他的形式可供选用，比如：

$\frac{dD}{dN}=\alpha {\left(\frac{{\sigma }_a}{1-{\mathrm{n\sigma }}_m}\right)}^m{(1-D)}^{-\beta }$

或者

$\frac{dD}{dN}=\frac{\beta }{E^m}\frac{{({\sigma }_{max}-{\sigma }_{th})}^m}{{(1-D)}^{(3/2)m}}$

其中，α、n、β、m为损伤参数，σ_a代表应力幅，σ_m代表平均应力。本文选用单轴加载的损伤演化方程进行实施例说明，对其他形式的演化方程本方法以下步骤依然适用不再详细说明。

(2)将损伤演化方程耦合传统有限元方法，得到给定损伤场时结构响应分析的损伤力学有限元分析列式：

$\left(\underset{e}{\Sigma }{\left[A_e\right]}^T\right[K_e\left]\right[A_e]-\underset{e}{\Sigma }{D}_e{\left[A_e\right]}^T[K_e\left]\right[A_e\left]\right)\left\{\delta \right\}=\left\{f\right\}$

其中，[A_e]为位移协调矩阵，[K_e]为单元刚度矩阵，D_e为单元损伤度，{δ}为总体位移列阵，{f}为载荷列向量；

有限元分析结构含损伤的响应时，需要使用耦合损伤度的损伤力学有限元方程，由上式可以看出，在已知单元的损伤度时，对于结构的位移响应分析只需要将单元刚度矩阵进行按损伤度折减后在进行常规有限元计算即可得到位移应力与应变，因此在有限元软件中通过每一步计算损伤度后给单元重新赋予折减后的刚度很容易实现；

(3)查询疲劳试验手册获取标准疲劳试验件的中值疲劳寿命数据拟合损伤演化方程中的损伤参量，将损伤演化方程在0到1上积分，对应一个单元破坏则是标准的S-N曲线的寿命值，然后用最小二乘法拟合β、α和p，整个拟和参数的详细流程如图2所示；

首先，通过试验和理论中值疲劳曲线确定参数p值与K_m值。理论中值疲劳曲线可采用下式来表示：

$({\sigma }_{\max }^{2p+2}-{\sigma }_{thm}^{2p+2})N_{fm}=K_m$

其中：

$K_m=\frac{p+1}{2p+3}\frac{{\left(2E\right)}^{p+1}{(1-D_{0m})}^{2p+3}}{\alpha }$

σ_thm＝σ_th0(1-D_0m)^0.5+β

式中，p、β为损伤参数，σ_max为外载荷最大应力，N_fm为理论疲劳寿命，σ_th0为无初始损伤情况对应的应力门槛值，σ_thm为中值疲劳曲线对应的应力门槛值，D_0m为中值疲劳曲线对应的初始损伤大小。

将理论中值疲劳曲线两边取对数写成如下形式