CN106529017A

CN106529017A - 一种高铬钢构件高温蠕变变形预测与蠕变损伤分析方法

Info

Publication number: CN106529017A
Application number: CN201610978437.1A
Authority: CN
Inventors: 王炯; 蔡晓丹
Original assignee: South China University of Technology SCUT
Current assignee: South China University of Technology SCUT
Priority date: 2016-11-08
Filing date: 2016-11-08
Publication date: 2017-03-22

Abstract

本发明公开了一种高铬钢构件高温蠕变变形预测与蠕变损伤分析方法，包括理论模型的建立、材料参数的确定、数值积分算法设计、有限元软件的二次开发、构件的蠕变损伤分析等。本发明方法相对于现有的技术方案进行了较大的改进，可以实现对高铬钢构件的蠕变行为和蠕变损伤的精确预测，因此，在超临界发电机组中高温高压构件的安全设计与剩余寿命评估领域将具有重要的应用价值。

Description

一种高铬钢构件高温蠕变变形预测与蠕变损伤分析方法

技术领域

本发明涉及超临界发电机组中高温高压构件的安全设计与剩余寿命评估的技术领域，尤其是指一种高铬钢构件高温蠕变变形预测与蠕变损伤分析方法。

背景技术

高铬钢(如P91钢、P92钢、E911钢等)材料作为一类新型抗蠕变合金材料，在高温高压条件下具有高蠕变强度、高硬度、耐腐蚀、易焊接等优良性能，故已在世界范围内被广泛应用于超临界电厂的建设中，并被视为制造下一代超超临界发电机组的理想材料。

在高温和长期加载条件下，高铬钢构件所面临的主要的运行风险之一是蠕变损伤和断裂。相关领域已对高铬钢材料的蠕变性能进行了系统的实验研究[1-3]。实验数据表明，高铬钢材料的蠕变性能具有下述主要特征：1)整个蠕变过程可分为三个阶段，其中第一阶段具有较高的蠕变速率并积累了较大的蠕变变形，需要特别关注[4]；2)在第一蠕变阶段中，蠕变速率与加载时间在Log-Log图中具有线性关系[5]；3)材料的最小蠕变速率与平均蠕变断裂时间满足线性Monkman-Grant关系[6]；4)高铬钢焊接件的蠕变强度远小于其原材料的蠕变强度，因此高铬钢焊接件的第Ⅳ型蠕变断裂对于发电机组设备的总体安全具有非常大的危害性[7]。

除实验研究之外，还需要建立合适的理论模型，对高铬钢构件的蠕变损伤进行模拟和预测，从而为构件的安全设计和剩余寿命预测提供理论依据。经典的蠕变法则(如Norton法则、Larson-Miller参数等)通常只关注于最小蠕变速率与应力、时间和温度的关系，无法满足构件安全设计的需求。近年来，科研人员也已提出了多个先进的本构模型对高铬钢材料的蠕变行为进行模拟。Hayhurst et al.[8]所建立的连续损伤模型中包含三个内部变量，分别表示材料的应变硬化、析出碳化物的体积和晶包间空隙损伤，可以对材料的蠕变过程进行比较全面的描述。Bauer et al.[9]为E911焊接管提出了一个修正的Graham-Walles蠕变法则，其中引入了一个损伤变量用于描述第三蠕变阶段。该模型可用于分析在内部压力和轴向加载条件下焊接管内部的等价应力和应变的分布。Hyde et al.[10]研究了两种蠕变连续损伤模型，即Kachanov模型和Liu-Murakami模型，通过比较后发现Liu-Murakami模型在模拟高铬钢材料的稳定性方面具有更好的稳定性。Naumenko et al.[11]在模型中引入了蠕变硬化变量和损伤变量，从而对蠕变的第一阶段和第三阶段进行了更好的描述。Chang et al.[12]基于Norton-Bailey时间硬化法则和Kachanov–Robotnov连续损伤理论建立了本构模型，用于对高铬钢材料的三个蠕变阶段进行描述，并考虑了多轴应力状态的影响。其他的本构模型可参见文献[13-15]。

为了满足高铬钢构件安全设计和剩余寿命预测的需求，上述模型还需要得到进一步改进。现有模型存在的主要问题包括：

1、模型的预测精度还需要得到进一步提高。特别是现有模型很少能够精确预测高铬钢材料在第一蠕变阶段蠕变速率的变化。由于第一蠕变阶段持续的时间较短，在本构建模过程中经常被忽略。而事实上，第一阶段的蠕变变形已构成了安全设计中构件的许可变形的主要部分[3]。另外，构件的局部区域的变形对于构件内部应力的分布具有很大的影响，从而会直接影响到构件的剩余使用寿命。所以，对于第一阶段蠕变变形的精确预测对于高铬钢构件的安全设计具有重要的意义。

2、对于较为复杂的本构模型，模型中材料参数的确定通常是一个难题。特别是对含有多个内部变量的本构模型，模型中不同因素的交互影响使得材料参数的确定变得非常困难。通常需要设计适当的方案确定参数的初值，再采用优化程序对参数值进行调整，以提高模型的预测精度[16]。对于不同的温度和应力加载条件，这一过程通常需要多次重复，从而消耗大量的计算资源。

除此之外，还需设计有效的数值积分算法对模型的本构发展方程组进行求解。材料模型还应被植入有限元程序(软件)中进行计算，从而实现对真实加载条件下高铬钢构件整体热力学行为的模拟和预测。

参考文献：

[1]European Creep Collaborative Committee,ECCC data sheets,OMMI 2005.

[2]National Research Institute for Metals(NIMS)Japan,Data sheets onthe elevated-temperature properties of 9Cr-1Mo-V-Nb steel tubes for boilersand heat exchanges(ASME SA-213/SA-213M Grade T91)and 9Cr-1Mo-V-Nb steelplates for boilers and pressure vessels(ASME SA-387/SA-387M Grade 91),1996.

[3]National Research Institute for Metals(NIMS)Japan,Data sheets onthe elevated-temperature properties of 9Cr-0.5Mo-1.8W-V-Nb steel tube forpower boilers(ASME SA-213/SA-213M Grade T92)and 9Cr-0.5Mo-1.8W-V-Nb steelpipe for high temperature service(ASME SA-335/SA-335M Grade P92),2012.

[4]L.Esposito,N.Bonora,Primary creep modeling based on the dependenceof the activation energy on the internal stress,Journal of Pressure VesselTechnology134(2012)061401.

[5]K.Naumenko,H.Altenbach,A.Kutschke,A combined model for hardening,softening,and damage processes in advanced heat resistant steels at elevatedtemperature,International Journal of Damage Mechanics 20(2011)578-597.

[6]K.Kimura,O.Kanemaru,K.Sawada,H.Kushima,Creep deformation andrupture strength property of ASME grades T/P92 Steels,38th MPA-Seminar,Oct.1-2,2012,Stuttgart.

[7]K.Laha,K.S.Chanderavathi,P.Parameswaranp,K.Bhanu Sankara Rao,TypeIV cracking susceptibility in weld joints of different grades of Cr-Moferritic steel,Metallurgical and Materials Transactions A,2009,40:386-397.

[8]R.J.Hayhurst,R.Mustata,D.R.Hayhurst,Creep constitutive equationsfor parent,Type IV,R-HAZ,CG-HAZ and weld material in the range 565-640℃forCr-Mo-V weldments,International Journal of Pressure Vessels and Piping 82(2005)137-144.

[9]M.Bauer,A.Klenk,K.Maile,E.Roos,Numerical investigations onoptimisation of weld creep performance in martensitic steels,In:Cerjak H,Bhadeshia HKDH,Kozeschnik E,editors.Mathematical Modelling of Weld Phenomena8(2007)409-423.

[10]T.H.Hyde,M.Saber,W.Sun,Testing and modelling of creep crackgrowth in compact tension specimens from a P91 weld at 650℃,EngineeringFracture Mechanics 77(2010)2946-2957.

[11]K.Naumenko,H.Altenbach,A.Kutschke,A combined model for hardening,softening,and damage processes in advanced heat resistant steels at elevatedtemperature,International Journal of Damage Mechanics 20(2011)578-597.

[12]Y.Chang,H.Xu,Y.Z.Ni,X.Lan,H.Y.Li,The effect of multiaxial stressstate on creep behavior and fracture mechanism of P92 steel,Material Scienceand Engineering A 636(2015)70-76.

[13]G.Lewis,K.M.Shaw,Creep constitutive model and component lifetimeestimation:the case of Niobium-modified 9Cr-1Mo steel weldments,Journal ofMaterials Engineering and Performance 20(2011)1310-1314.

[14]T.Ogata,T.Sakai,M.Yaguchi,Damage characterization of a P91 steelweldment under uniaxial and multiaxial creep,Materials Science andEngineering A510-511(2009)238-243.

[15]M.Basirat,T.Shrestha,G.Potirniche,I.Charit,K.Rink,A study of thecreep behavior of modified 9cr-1mo steel using continuum-damage modeling,International Journal of Plasticity 37(2012)95-107.

[16]Y.P.Gong,C.J.Hyde,W.Sun,T.H.Hyde,Determination of materialproperties in the Chaboche unified viscoplasticity model,Proceedings of theInstitution of Mechanical Engineers,Part L:Journal of Materials:Design andApplications 224(2010)19-29.

发明内容

本发明的目的在于克服现有技术的缺点和不足，提供一种高铬钢构件高温蠕变变形预测与蠕变损伤分析方法，该方法可以实现对高铬钢构件的蠕变行为和蠕变损伤的精确预测，因此，在超临界发电机组中高温高压构件的安全设计与剩余寿命评估领域将具有重要的应用价值。

为实现上述目的，本发明所提供的技术方案为：一种高铬钢构件高温蠕变变形预测与蠕变损伤分析方法，包括以下步骤：

1)建立新的材料本构模型，即理论模型，对高铬钢材料在整个蠕变过程中的热力学响应进行模拟，并对材料的蠕变损伤进行分析；其中，建立的理论模型如下：

总应变张量：

胡克定律：σ＝σ^vol+S,S＝2GE^e

总蠕变速率：m_H＝3S/(2||σ||_H)

第一蠕变阶段的蠕变速率：

第二蠕变阶段的蠕变速率：

P_s＝a₀+a₁log(||σ||_H)+a₂log(||σ||_H)²+a₃log(||σ||_H)³

第三蠕变阶段的蠕变速率：

P_t＝A₀+A₁log(||σ||_H)+A₂log(||σ||_H)²+A₃log(||σ||_H)³

上述模型中，ε为总应变张量，为弹性应变张量ε^e和蠕变应变张量ε^c之和，ε^e进一步分解为体积应变部分和偏应变部分E^e之和，其中Tr(·)表示张量的迹，I为单位张量；σ为应力张量，它与弹性应变张量ε^e满足胡克定律，相应于ε^e的分解，σ也能够分解为体积应力部分σ^vol和偏应力部分S，其中K和G分别为体积弹性模量和剪切弹性模量；总蠕变速率记为其方向为m_H，总蠕变速率的值等于第一阶段蠕变速率因子e^pc、第二阶段的最小蠕变速率e^sc和第三阶段蠕变速率因子e_tc的乘积，其中||σ||_H为von Mises等效应力；为了模拟第一蠕变阶段的材料变形硬化效应，蠕变速率因子e^pc的定义依赖于硬化变量H，硬化变量的变化率与第一蠕变阶段的持续时间t_h成反比，而t_h能够进一步表示为平均蠕变断裂时间t_r的函数，其中A_pc,n_pc,C_h,n_h是相关的材料参数；第二蠕变阶段的最小蠕变速率e^sc通过Larson-Miller参数P_s进行定义，而P_s能够进一步写为von Mises等效应力||σ||_H的函数，其中T为绝对温度值，a₀,a₁,a₂,a₃,C_s为相关的材料参数；第三蠕变阶段的蠕变速率因子e^tc基于损伤变量D^c进行定义，损伤变量的变化率依赖于平均蠕变断裂时间t_r，而t_r的定义依赖于另一个Larson-Miller参数P_t，其中A₀,A₁,A₂,A₃,C_t,n_tc为相关的材料参数；

2)基于加载的实验数据，提出模型中材料参数的确定方案，实现对不同温度和应力加载条件下高铬钢材料热力学响应的精确预测；其中，材料参数的确定如下：

首先，根据高铬钢材料在单轴拉伸实验的应力-应变曲线的斜率，即线性变形阶段，确定材料的杨氏模量E；对于金属材料，泊松比取为ν＝0.3；由此，弹性模量K和G的值根据下式确定

为了确定与e^sc相关的材料参数a₀,a₁,a₂,a₃,C_s，考虑高铬钢材料在单轴加载蠕变实验中所测量的不同温度和应力条件下的最小蠕变速率值；基于实验结果，并采用最小二乘拟合算法，即可确定材料参数a₀,a₁,a₂,a₃,C_s的值；需要指出的是a₀,a₁,a₂,a₃,C_s的值不依赖于温度和应力的大小；

为了确定与e^tc相关的材料参数A₀,A₁,A₂,A₃,C_t，考虑高铬钢材料在单轴加载蠕变实验中所测量的不同温度和应力条件下的平均蠕变断裂时间；基于实验结果，并采用最小二乘拟合算法，即可确定材料参数A₀,A₁,A₂,A₃,C_t的值；需要指出的是A₀,A₁,A₂,A₃,C_t的值同样不依赖于温度和应力的大小；参数n_tc的值能够进行调整，从而使得第三阶段的蠕变曲线得到更好的模拟；

为了确定与e^pc相关的材料参数C_h和n_h的值，考虑高铬钢材料在单轴加载蠕变实验中第一蠕变阶段时间t_h和平均蠕变断裂时间t_r的关系；根据实验数据，在log(t_r)-log(t_h)图表中，不同应力条件下的(log(t_r),log(t_h))点落在同一条直线上；因此，理论模型中假设t_h和t_r满足幂函数的关系式，其中的材料参数C_h和n_h能够通过与实验数据的拟合得以确定；

为了确定与e^pc相关的材料参数A_pc和n_pc的值，考虑高铬钢材料在单轴加载蠕变实验中第一蠕变阶段的应变速率和最小蠕变速率e^sc的关系；推导得到：

另一方面，根据实验数据，在图表中，不同时刻和不同应力条件下的数据点落在同一条直线上，这也证明了关系式(2)的正确性；基于实验数据对关系式(2)进行拟合，即可确定参数A_pc和n_pc的值；

3)设计数值积分算法，并进行有限元软件的二次开发，对各种高铬钢构件在不同加载条件下的热力学行为进行模拟；其中，

数值积分算法的设计如下：

采用后向Euler数值积分的方法对理论模型的本构发展方程组进行求解，对于应变控制的问题，设计算法如下：

3.1.1)在时刻t_n已知状态变量在当前时刻t＝t_n+Δt，已知总应变ε；

3.1.2)通过胡克定律计算体积应力张量σ^vol，并推导蠕变应变张量ε^c和偏应力张量S的关系式设定当前时刻状态变量的初值X＝{S,H,D}；

3.1.3)构建剩余方程组Π＝{R_S,R_H,R_D}^T＝0^T，其中

3.1.4)收敛性检验：如果计算收敛，跳出程序；

3.1.5)计算切向算子其中

3.1.6)求解线性方程组更新状态变量X+ΔX→X，回到步骤3.1.3)；

有限元软件的二次开发如下：

为了模拟超临界发电机组中的高铬钢构件在真实加载条件下的热力学行为并分析其蠕变损伤，需要将上述材料模型植入有限元软件中，开展有限元模拟；选取软件ANSYS进行有限元计算，基于ANSYS Usermat子程序完成了二次开发，实现了对高铬钢构件整体热力学行为的模拟；具体的实施步骤如下：

3.2.1)基于ANSYS Usermat子程序的框架结构编写C语言程序，实现材料模型的数值积分算法；

3.2.2)将编写的C语言程序在ANSYS custom平台进行编译，生成用户定义的ANSYS可执行文件；

3.2.3)针对不同的问题，编写ANSYS APDL文件，调用用户定义的ANSYS可执行文件，进行有限元计算；

4)根据有限元模拟的结果对高铬钢构件的蠕变损伤进行分析，为构件的安全设计和剩余寿命评估提供理论依据。

本发明与现有技术相比，具有如下优点与有益效果：

1、完善的理论模型

本发明所建立的理论模型全面考虑了高铬钢材料在蠕变过程中的热力学行为特征，包括蠕变第一阶段的材料硬化效应、第二阶段的最小蠕变速率和第三阶段的蠕变损伤断裂。另外，由于该模型在三维框架下建立，故其可用于描述高铬钢材料在多轴加载条件下的热力学响应，从而为模拟真实高铬钢构件在复杂加载条件下的热力学行为奠定了基础。

2、成熟的材料参数确定方案

基于理论模型的特征，本发明提出了完整的材料参数确定方案。该方案只需采用高铬钢材料单轴蠕变实验的基本数据(如最小蠕变速率、平均蠕变断裂时间、第一阶段蠕变速率变化曲线等)，这些数据可在ECCC和NIMS的数据表中查询得到。由该方案所确定的材料参数无需再做优化便可直接用于高铬钢材料力学响应的预测。

3、对高铬钢材料局部热力学响应的精确预测

本发明所提出的技术方案(包括理论模型、参数确定方案、数值积分算法设计等)可以实现对高铬钢材料局部热力学响应的精确预测。以ASME Grade 92钢为例，图4a‐图7b展示了本方案的预测结果与实验结果的比较。

4、对高铬钢构件整体热力学行为及其蠕变损伤的有效模拟和预测

基于有限元软件的二次开发，本发明可以实现对高铬钢构件在复杂加载条件下的整体热力学行为的模拟，并对构件中关键部位(如应力最大区域或应变最大区域)的蠕变损伤进行分析。

为了举例说明上述效果，考虑如图8a所示高铬钢集管构件。由于构件形状和加载方式的对称性，有限元模拟中只需考虑如图8b所示的构件部分。考虑该集管构件承受内部蒸汽压力的作用。在构件的内表面(图8b所示P面)施加均匀压力，在构件管道的截面(图8b所示S面)上设置对称边界条件。另外，在喷嘴管道上施加轴向拉力以保证其平衡。整个加载过程分为两个阶段，第一阶段持续时间为1小时，压强由0MPa增加到15MPa。第二阶段持续时间为10000小时，压强保持在15MPa。

通过调用用户开发的ANSYS Usermat子程序，可以对上述高格钢构件的热力学行为进行有限元模拟。其中，材料参数和加载方式均在ANSYS APDL文件中进行设定。根据有限元计算的结果，可以确定构件中应力、应变等状态变量的分布(图9a‐图10b)。

附图说明

图1为本发明的技术路线图。

图2为ASME Grade 92钢在650℃和不同应力条件下的单轴蠕变实验中的logt_r-logt_h关系图及材料参数C_h和n_h的拟合结果。

图3为T92钢在650℃和不同应力条件下的单轴蠕变实验中的关系图及材料参数A_pc和n_pc的拟合结果。

图4a为ASME Grade 92钢的最小蠕变速率的预测结果(实线)与实验结果(点线)的比较(实验结果取自文献[6])。

图4b为ASME Grade 92钢的平均蠕变断裂时间的预测结果(实线)与实验结果(点线)的比较(实验结果取自文献[6])。

图5为ASME Grade 92钢在蠕变第一阶段蠕变速率变化曲线的预测结果(实线)与实验结果(点线)的比较(实验结果取自文献[6])。

图6为ASME Grade 92钢在整个蠕变过程中的蠕变速率变化曲线的预测结果与实验结果的比较(实验结果取自文献[6])。

图7a为ASME Grade 92钢在650℃的蠕变曲线的预测结果(黑线)与实验结果(点)的比较(实验结果取自文献[6])。

图7b为ASME Grade 92钢在700℃的蠕变曲线的预测结果(黑线)与实验结果(点)的比较(实验结果取自文献[6])。

图8a为计算示例选取的高铬钢集管构件几何尺寸(单位：mm)及其所承受加载方式图之一。

图8b为计算示例选取的高铬钢集管构件几何尺寸(单位：mm)及其所承受加载方式图之二。

图9a为在第一加载阶段末高铬钢构件中Von Mises等效应力的分布。

图9b为在第二加载阶段末高铬钢构件中Von Mises等效应力的分布。

图10a为在第一加载阶段末高铬钢构件中等效非弹性应变的分布。

图10b为在第二加载阶段末高铬钢构件中等效非弹性应变的分布。

图11为高铬钢构件的关键材料点处的蠕变损伤变量在整个蠕变过程中的演化曲线。

具体实施方式

下面结合具体实施例对本发明做进一步的说明。

本发明旨在为高铬钢构件高温蠕变变形预测与蠕变损伤分析提出系统的解决方案。方案的内容包括理论模型的建立、材料参数的确定、数值积分算法设计、有限元软件的二次开发、构件的蠕变损伤分析等。该方案将通过图1所示的技术路线图实施，具体包括以下步骤：

1)建立新的材料本构模型，即理论模型，对高铬钢材料在整个蠕变过程(包含三个阶段)中的热力学响应进行模拟，并对材料的蠕变损伤进行分析。

建立的理论模型如下：

总应变张量：

胡克定律：σ＝σ^vol+S,S＝2GE^e

总蠕变速率：m_H＝3S/(2||σ||_H)

第一蠕变阶段的蠕变速率：

第二蠕变阶段的蠕变速率：

P_s＝a₀+a₁log(||σ||_H)+a₂log(||σ||_H)²+a₃log(||σ||_H)³

第三蠕变阶段的蠕变速率：

P_t＝A₀+A₁log(||σ||_H)+A₂log(||σ||_H)²+A₃log(||σ||_H)³

上述模型中，ε为总应变张量，为弹性应变张量ε^e和蠕变应变张量ε^c之和，ε^e可进一步分解为体积应变部分和偏应变部分E^e之和，其中Tr(·)表示张量的迹，I为单位张量；σ为应力张量，它与弹性应变张量ε^e满足胡克定律，相应于ε^e的分解，σ也能够分解为体积应力部分σ^vo^l和偏应力部分S，其中K和G分别为体积弹性模量和剪切弹性模量；总蠕变速率记为其方向为m_H，总蠕变速率的值等于第一阶段蠕变速率因子e^pc、第二阶段的最小蠕变速率e^sc和第三阶段蠕变速率因子e^tc的乘积，其中||σ||_H为von Mises等效应力；为了模拟第一蠕变阶段的材料变形硬化效应，蠕变速率因子e^pc的定义依赖于硬化变量H，硬化变量的变化率与第一蠕变阶段的持续时间t_h成反比，而t_h能够进一步表示为平均蠕变断裂时间t_r的函数，其中A_pc,n_pc,C_h,n_h是相关的材料参数；第二蠕变阶段的最小蠕变速率e^sc通过Larson-Miller参数P_s进行定义，而P_s能够进一步写为von Mises等效应力||σ||_H的函数，其中T为绝对温度值，a₀,a₁,a₂,a₃,C_s为相关的材料参数；第三蠕变阶段的蠕变速率因子e^tc基于损伤变量D_c进行定义，损伤变量的变化率依赖于平均蠕变断裂时间t_r，而t_r的定义依赖于另一个Larson-Miller参数P_t，其中A₀,A₁,A₂,A₃,C_t,n_tc为相关的材料参数。

2)基于简单加载的实验数据，提出模型中材料参数的确定方案，实现对不同温度和应力加载条件下高铬钢材料热力学响应的精确预测。

本发明所提出的理论模型包含下述17个材料参数：

①弹性模量：K,G

②与e^sc相关的材料参数：a₀,a₁,a₂,a₃,C_s

③与e^tc相关的材料参数：A₀,A₁,A₂,A₃,C_t,n_tc

④与e^pc相关的材料参数：A_pc,n_pc,C_h,n_h

为了确定材料参数的值，我们采用如下的材料参数确定方案：

首先，根据高铬钢材料在单轴拉伸实验的应力-应变曲线的斜率(线性变形阶段)，可以确定材料的杨氏模量E。对于金属材料，泊松比可取为ν＝0.3。由此，弹性模量K和G的值可以根据下式确定

为了确定与e^sc相关的材料参数a₀,a₁,a₂,a₃,C_s，考虑高铬钢材料在单轴加载蠕变实验中所测量的不同温度和应力条件下的最小蠕变速率值(见文献[2,3,6])。基于实验结果，并采用最小二乘拟合算法，即可确定材料参数a₀,a₁,a₂,a₃,C_s的值。需要指出的是a₀,a₁,a₂,a₃,C_s的值不依赖于温度和应力的大小。

为了确定与e^tc相关的材料参数A₀,A₁,A₂,A₃,C_t，考虑高铬钢材料在单轴加载蠕变实验中所测量的不同温度和应力条件下的平均蠕变断裂时间(见文献[1-3,6])。基于实验结果，并采用最小二乘拟合算法，即可确定材料参数A₀,A₁,A₂,A₃,C_t的值。需要指出的是A₀,A₁,A₂,A₃,C_t的值同样不依赖于温度和应力的大小。参数n_tc的值可以进行适当的调整，从而使得第三阶段的蠕变曲线可以得到更好的模拟(通常可以取n_tc＝1)。

为了确定与e^pc相关的材料参数C_h和n_h的值，考虑高铬钢材料在单轴加载蠕变实验中第一蠕变阶段时间t_h和平均蠕变断裂时间t_r的关系。根据实验数据(见文献[2,3,6])可以发现，在log(t_r)-log(t_h)图表中，不同应力条件下的(log(t_r),log(t_h))点落在同一条直线上。因此，理论模型中假设t_h和t_r满足幂函数的关系式(见表1第四项)，其中的材料参数C_h和n_h可以通过与实验数据的拟合得以确定(见图2示例)。

为了确定与e^pc相关的材料参数A_pc和n_pc的值，考虑高铬钢材料在单轴加载蠕变实验中第一蠕变阶段的应变速率和最小蠕变速率e^sc的关系。由表1第4项的关系式出发，通过简单的推导可得

另一方面，由实验数据可知在图表中，不同时刻和不同应力条件下的数据点落在同一条直线上(见图3)，这也证明了关系式(2)的正确性。基于实验数据对关系式(2)进行拟合，即可确定参数A_pc和n_pc的值(见图3示例)。

综上所述，尽管理论模型包含的材料参数较多。但这些材料参数可以通过完善的方案得以确定，且参数的确定仅需提供高铬钢材料在单轴加载蠕变试验中的数据。后续计算可以表明，采用上述方案所确定的材料参数可以对高铬钢材料的蠕变响应进行精确的预测。

3)设计数值积分算法，并进行有限元软件的二次开发，对各种高铬钢构件在不同加载条件下的热力学行为进行模拟。

本发明专利采用后向Euler数值积分的方法对理论模型的本构发展方程组进行求解。对于应变控制的问题，设计算法如下：

3.1.1)在时刻t_n已知状态变量在当前时刻t＝t_n+Δt，已知总应变ε。

3.1.2)通过胡克定律计算体积应力张量σ^vol，并推导蠕变应变张量ε^c和偏应力张量S的关系式设定当前时刻状态变量的初值X＝{S,H,D}。

3.1.3)构建剩余方程组Π＝{R_S,R_H,R_D}^T＝0^T，其中

3.1.4)收敛性检验：如果计算收敛，跳出程序。

3.1.5)计算切向算子其中

3.1.6)求解线性方程组更新状态变量X+ΔX→X，回到步骤3.1.3)。

为了模拟超临界发电机组中的高铬钢构件在真实加载条件下的热力学行为并分析其蠕变损伤，需要将上述材料模型植入有限元程序(软件)中，开展有限元模拟。这一任务可以采用不同的有限元程序(软件)来实现。

在此，我们选取商业软件ANSYS进行有限元计算，基于ANSYS Usermat子程序完成了二次开发，实现了对高铬钢构件整体热力学行为的模拟。具体的实施步骤如下：

3.2.1)基于ANSYS Usermat子程序的框架结构编写C语言程序，实现材料模型的数值积分算法。

3.2.2)将编写的C语言程序在ANSYS custom平台进行编译，生成用户定义的ANSYS可执行文件。

3.2.3)针对不同的问题，编写ANSYS APDL文件，调用用户定义的ANSYS可执行文件，进行有限元计算。

4)根据有限元模拟的结果对高铬钢构件的蠕变损伤进行分析，为构件的安全设计和剩余寿命评估提供理论依据。例如，选取第一加载阶段末具有最大von Mises等效应力的材料点，该材料点处的蠕变损伤变量在整个加载过程中的演化曲线如图11所示。根据该曲线可以预测不同时刻该点处的蠕变损伤的程度，并判断其的发生蠕变断裂的时间。

以上所述实施例只为本发明之较佳实施例，并非以此限制本发明的实施范围，故凡依本发明之形状、原理所作的变化，均应涵盖在本发明的保护范围内。

Claims

1.一种高铬钢构件高温蠕变变形预测与蠕变损伤分析方法，其特征在于：包括以下步骤：

总应变张量：

胡克定律：σ＝σ^vol+S,S＝2GE^e

总蠕变速率：m_H＝3S/(2||σ||_H)

第一蠕变阶段的蠕变速率：

e^{p c} = 1 + A_{p c} < \frac{1}{H} - 1 >^{n_{p c}}, \overset{\cdot}{H} = 1 / t_{h}, t_{h} = C_{h} {t_{r}}^{n_{h}}

第二蠕变阶段的蠕变速率：

e^{s c} (| | σ | |_{H}, T) = 10^{P_{s} / T - C_{s /}} 3600,

P_s＝a₀+a₁ log(||σ||_H)+a₂ log(||σ||_H)²+a₃ log(||σ||_H)³

第三蠕变阶段的蠕变速率：

e^{t c} (D^{c}) = {(1 - D^{c})}^{- n_{t c}}, {\overset{\cdot}{D}}^{c} = 1 / t_{r},

t_{r} (| | σ | |_{H}, T) = 3600 \times 10^{P_{t} / T - C_{t}},

P_t＝A₀+A₁ log(||σ||_H)+A₂ log(||σ||_H)²+A₃ log(||σ||_H)³

上述模型中，ε为总应变张量，为弹性应变张量ε^e和蠕变应变张量ε^c之和，ε^e进一步分解为体积应变部分和偏应变部分E^e之和，其中Tr(·)表示张量的迹，I为单位张量；σ为应力张量，它与弹性应变张量ε^e满足胡克定律，相应于ε^e的分解，σ也能够分解为体积应力部分σ^vol和偏应力部分S，其中K和G分别为体积弹性模量和剪切弹性模量；总蠕变速率记为其方向为m_H，总蠕变速率的值等于第一阶段蠕变速率因子e^pc、第二阶段的最小蠕变速率e^sc和第三阶段蠕变速率因子e^tc的乘积，其中||σ||_H为von Mises等效应力；为了模拟第一蠕变阶段的材料变形硬化效应，蠕变速率因子e^pc的定义依赖于硬化变量H，硬化变量的变化率与第一蠕变阶段的持续时间t_h成反比，而t_h能够进一步表示为平均蠕变断裂时间t_r的函数，其中A_pc,n_pc,C_h,n_h是相关的材料参数；第二蠕变阶段的最小蠕变速率e^sc通过Larson-Miller参数P_s进行定义，而P_s能够进一步写为vonMises等效应力||σ||_H的函数，其中T为绝对温度值，a₀,a₁,a₂,a₃,C_s为相关的材料参数；第三蠕变阶段的蠕变速率因子e^tc基于损伤变量D^c进行定义，损伤变量的变化率依赖于平均蠕变断裂时间t_r，而t_r的定义依赖于另一个Larson-Miller参数P_t，其中A₀,A₁,A₂,A₃,C_t,n_tc为相关的材料参数；

K = \frac{E}{3 (1 - 2 v)}, G = \frac{E}{2 (1 + v)} . - - - (1)

\frac{| | {\overset{\cdot}{ϵ}}^{c} | |_{H}}{e^{s c}} - 1 = A_{p c} < \frac{t_{h}}{t} - 1 >^{n_{p c}}, 0 < t < t_{h} . - - - (2)

数值积分算法的设计如下：

3.1.3)构建剩余方程组Π＝{R_S,R_H,R_D}^T＝0^T，其中

R_{S} = - \frac{S - S_{n}}{2 G} + E - E_{n} - e^{p c} (H) \cdot e^{s c} (| | σ | |_{H}) \cdot e^{t c} (D^{c}) \cdot m_{H} Δ t

R_{H} = H - H_{n} - \frac{Δ t}{t_{h} (| | σ | |_{H})}

R_{D} = D^{c} - D_{n}^{c} - \frac{Δ t}{t_{r} (| | σ | |_{H})}

3.1.4)收敛性检验：如果计算收敛，跳出程序；

3.1.5)计算切向算子其中

A = {[\begin{matrix} \frac{\partial R_{S}}{\partial S} & \frac{\partial R_{S}}{\partial H} & \frac{\partial R_{S}}{\partial D} \\ \frac{\partial R_{H}}{\partial S} & \frac{\partial R_{H}}{\partial H} & \frac{\partial R_{H}}{\partial D} \\ \frac{\partial R_{D}}{\partial S} & \frac{\partial R_{D}}{\partial H} & \frac{\partial R_{D}}{\partial D} \end{matrix}]}_{8 \times 8}

3.1.6)求解线性方程组更新状态变量X+ΔX→X，回到步骤3.1.3)；

有限元软件的二次开发如下：