CN104809273A - 一种描述蠕变变形的方法 - Google Patents

Abstract

一种描述蠕变变形的方法，它是一种能完整描述3个阶段蠕变变形的方法，步骤如下：一、将试验得到的不同温度和应力下的蠕变曲线归一化；二、对曲线的第2阶段进行拟合，得到η₁，η₂，η₃；三、拟合得到η₄，η₅；四、将η_i该i＝1,2,3,4,5，表示为温度和应力的函数，将η_i该i＝1,2,3,4,5表示为无量纲应力σ/σ_0.2与无量纲温度T/T_m的函数；五、将该方法与有限元软件结合，编写usercreep子程序；六、选择适当的适用于变载情况下的蠕变模型；七、建立实际结构的有限元模型，进行蠕变变形计算和应力松弛行为分析；所发明的蠕变变形方法，解决了当前常用的蠕变模型描述蠕变曲线的不足，同时针对变载情况，在子程序中实现3种变载模型。

Description

一种描述蠕变变形的方法

一.技术领域

本发明提供一种能完整描述3个阶段蠕变变形的模型及将其应用于实际结构的蠕变计算的方法，属于高温结构蠕变模拟技术领域。

二.背景技术

蠕变是指在温度、载荷不变的条件下，材料的变形也会随时间增长而逐渐增大的现象，而且这种变形即使在应力小于屈服极限时仍具有不可逆的变形性质。典型的蠕变变形可以分为3个阶段：第1阶段由于变形引起加工硬化导致蠕变速率随时间不断降低，称为初始蠕变阶段；第2阶段为直线，这是由于加工硬化与回复软化过程达到动态平衡而导致蠕变速率保持不变，即为稳态蠕变阶段；第3阶段蠕变速率随时间增大直至断裂，称为加速蠕变阶段。

当前工程中常用的蠕变模型有时间硬化模型和应***化模型等，但这类模型只能模拟蠕变第1阶段或者前两个阶段(有限元软件ANSYS和ABAQUS提供的蠕变模型均不能模拟蠕变第3阶段)；1985年Evans和Wilshire在其著作中提出一种能模拟完整的蠕变变形的模型—θ投射法，虽然θ投射法可以模拟蠕变变形的3个阶段，但它并不注重描述稳态蠕变阶段，即它所描述的整个蠕变变形中，并不存在蠕变速率为常数的阶段；一些粘塑性本构模型虽能较好地描述循环载荷下的塑形变形，而在描述蠕变第3阶段方面必须引入损伤参量对模型进行修正。这些模型用于试样在试验载荷下变形的描述较精准，但由于计算较为复杂，目前还不易用于实际结构(如涡轮叶片和轮盘)的蠕变分析。

因此，仍有必要发展能够完整描述蠕变曲线的3个阶段的方法，同时能与有限元结合用于实际结构的蠕变变形计算及应力松弛问题的分析。

三.发明内容

1.发明目的

本发明提供一种能完整描述3个阶段蠕变变形的模型及将其应用于实际结构的蠕变计算的方法，达到能计算实际结构3个阶段的蠕变变形和模拟应力松弛行为的目的，解决现有技术的不足。

2.技术方案

本发明提供一种描述蠕变变形的方法，它能完整描述3个阶段(第1阶段由于变形引起加工硬化导致蠕变速率随时间不断降低，称为初始蠕变阶段；第2阶段蠕变速率不随时间变化，即稳态蠕变阶段；第3阶段蠕变速率随时间增大直至断裂，称为加速蠕变阶段)蠕变变形的方法，同时通过将该方法与有限元软件ANSYS结合，编写usercreep子程序，用于实际结构的蠕变变形计算。同时提出3种适用于变载情况的模型，用于计算应力松弛行为。本发明用归一化参数描述蠕变变形，将蠕变应变表示为：ε_c＝η₁ζ^η4+η₂ζ+η₃ζ^η5或分别称之为表达式1和表达式2，其中3项分别表示蠕变的3个阶段。表达式中各参数含义如下：ζ＝t/t_c为无量纲时间，t_c为给定温度和应力下的持久寿命，则ζ∈[0,1]；η_i(i＝1,2,3,4,5)为材料参数，η₁，η₂，η₃，分别为蠕变3个阶段的蠕变量，η₄，η₅分别控制蠕变第1阶段和第3阶段的变化快慢，且η₅>1。

本发明一种描述蠕变变形的方法，它是一种能完整描述3个阶段蠕变变形的方法，其具体步骤如下：

步骤一：将试验得到的不同温度和应力下的蠕变曲线归一化；

用试验得到的不同温度和应力下的蠕变曲线的时间坐标除以该温度和应力下的持久断裂时间，则所有曲线的横坐标均为归一化时间坐标ζ＝t/t_c，ζ∈[0,1]；

步骤二：对曲线的第2阶段进行拟合，得到η₁，η₂，η₃；

对每条曲线第2阶段进行拟合，若曲线第2阶段较为明显，则得到直线方程y＝kx+b，否则，根据最小蠕变率点得到直线方程；得到的一组k值则为不同应力和温度下的一组η₂值，即归一化坐标下的稳态蠕变应变率或最小蠕变率；拟合直线中的b值即η₁，继而可得到η₃＝ε_r-η₁-η₂，ε_r为断裂时的蠕变应变；

步骤三：拟合得到η₄，η₅；

然后由已得到的参数η₁，η₂，η₃，根据表达式或对每条曲线进行拟合(可采用matlab等软件进行最小二乘拟合)，得到η₄，η₅；

步骤四：将η_i(i＝1,2,3,4,5)表示为温度和应力的函数；

将η_i(i＝1,2,3,4,5)表示为无量纲应力σ/σ_0.2与无量纲温度T/T_m的函数；

步骤五：将该方法与有限元软件结合，编写usercreep子程序；

利用通用有限元软件ANSYS中usercreep子程序，将步骤一～步骤四得到的表达式写入子程序中，以达到计算实际结构蠕变变形的目的；对所编写的子程序进行编译连接后，在ANSYS主程序中调用子程序即可用所发明的方法对实际结构进行蠕变变形计算；

步骤六：选择适当的适用于变载情况下的蠕变模型；

在usercreep子程序中实现3种不同的适用于变载情况下蠕变模型：时间硬化模型，相对时间硬化模型，应***化模型；时间硬化模型为当载荷在t时刻变化(由σ₁,T₁变到σ₂,T₂)，t时刻后的蠕变曲线由σ₂,T₂状态下t时刻后的蠕变应变曲线上下平移得到；相对时间硬化模型为当载荷在t时刻变化(由σ₁,T₁变到σ₂,T₂)，t时刻后的蠕变曲线由σ₂,T₂状态下t₂＝t×t_c,1/t_c,2时刻后的蠕变曲线平移得到，t_c,1和t_c,2分别表示σ₁,T₁和σ₂,T₂状态下的持久寿命；应***化模型为载荷在t时刻变化(由σ₁,T₁变到σ₂,T₂)，t时刻后的蠕变曲线由σ₂,T₂状态下产生与前一状态相同的蠕变应变所对应的时间之后曲线左右平移得到；

步骤七：建立实际结构的有限元模型，进行蠕变变形计算和应力松弛行为分析；

对实际结构在CAE前处理软件中建立有限元模型；由于应力松弛行为属于变载条件下的蠕变行为，能通过所编写的usercreep子程序进行计算分析；在ANSYS主程序中调用所编写的usercreep子程序，同时输入所求得模型参数值，对实际结构进行有限元数值模拟，计算蠕变变形和应力松弛行为。

其中，在步骤二中所述的“对每条曲线第2阶段进行拟合，若曲线第2阶段较为明显，则得到直线方程y＝kx+b，否则，根据最小蠕变率点得到直线方程”，其方法如下：利用第二阶段的试验数据点，可采用excel或matlab等软件对该数据点进行一次多项式拟合；若蠕变曲线的第2阶段并不明显，可采用多次多项式对蠕变曲线进行拟合，然后对多次多项式求导，导数值最小的点即为最小蠕变率点，导数值即η₂。

其中，在步骤四中所述的“将η_i(i＝1,2,3,4,5)表示为温度和应力的函数”，具体做法如下：选取特定函数将无量纲应力σ/σ_0.2与无量纲温度T/T_m表示为η_i＝f(σ/σ_0.2,T/T_m)，例如：

${\mathrm{\η}}_i=a_i+b_i\frac{T}{T_m}+c_i\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}}+d_i\frac{T}{T_m}\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}}$ 或 $\ln {\mathrm{\η}}_i=a_i+b_i\frac{T}{T_m}+c_i\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}}+d_i\frac{T}{T_m}\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}},$

其中，a_i，b_i，c_i，d_i(i＝1,2,3,4,5)为材料相关系数，可由进行最小二乘拟合(可采用matlab等软件)得到；通过一组a_i，b_i，c_i，d_i的系数值，则可求得任意温度和应力下的η_i，继而得到该温度和应力下的蠕变曲线。

其中，在步骤五中所述的“usercreep子程序”，其编写所需的输出变量为：蠕变应变增量delcr、蠕变应变增量对等效应力的导数dcrda(1)，蠕变应变增量对蠕变应变的导数dcrda(2)；表达式2的各输出变量为：

蠕变应变增量delcr：

$\mathrm{delcr}={\stackrel{\·}{\mathrm{\ϵ}}}_c\mathrm{\Δt}=\frac{1}{t_c}({\mathrm{\η}}_1{\mathrm{\η}}_4{e}^{{-\mathrm{\η}}_4\mathrm{\ζ}}+{\mathrm{\η}}_2+{\mathrm{\η}}_3{\mathrm{\η}}_5{\mathrm{\ζ}}^{{\mathrm{\η}}_5-1})\mathrm{\Δt}$

蠕变应变增量对等效应力的导数dcrda(1)：

$\begin{array}{c}\mathrm{dcrda}\left(1\right)=\frac{\∂\left({\stackrel{\·}{\mathrm{\ϵ}}}_c\mathrm{\Δt}\right)}{\∂\mathrm{\σ}}\\ =\frac{1}{t_c}\frac{1}{{\mathrm{\σ}}_{0.2}}[{\mathrm{\η}}_1{\mathrm{\η}}_4{e}^{{-\mathrm{\η}}_4\mathrm{\ζ}}(q_1+q_4-q_4{\mathrm{\η}}_4\mathrm{\ζ})+q_2{\mathrm{\η}}_2+{\mathrm{\η}}_3{\mathrm{\η}}_5{\mathrm{\ζ}}^{{\mathrm{\η}}_5-1}(q_3+q_5+q_5{\mathrm{\η}}_5\ln \mathrm{\ζ})]\mathrm{\Δt}\end{array}$

式中， $q_i=c_i+d_i\frac{T}{T_m}(i=\mathrm{1,2,3,4,5});$

蠕变应变增量对蠕变应变的导数dcrda(2)：

$\mathrm{dcrda}\left(2\right)=\frac{\∂\left({\stackrel{\·}{\mathrm{\ϵ}}}_c\mathrm{\Δt}\right)}{\∂{\mathrm{\ϵ}}_c}=\frac{d{\stackrel{\·}{\mathrm{\ϵ}}}_c}{\mathrm{dt}}\left(\frac{1}{\frac{d{\mathrm{\ϵ}}_c}{\mathrm{dt}}}\right)\mathrm{\Δt}=\frac{{-\mathrm{\η}}_1{\mathrm{\η}}_4^2{e}^{{-\mathrm{\η}}_4\mathrm{\ζ}}+{\mathrm{\η}}_3{\mathrm{\η}}_5({\mathrm{\η}}_5-1){\mathrm{\ζ}}^{{\mathrm{\η}}_5-2}}{t_c({\mathrm{\η}}_1{\mathrm{\η}}_4{e}^{{-\mathrm{\η}}_4\mathrm{\ζ}}+{\mathrm{\η}}_2+{\mathrm{\η}}_3{\mathrm{\η}}_5{\mathrm{\ζ}}^{{\mathrm{\η}}_5-1})}\mathrm{\Δt}.$

其中，在步骤六中所述的“变载情况下的蠕变模型”，在子程序中实现方法为：

在子程序中，根据t₂＝t×t_c,1/t_c,2，可计算得到相对时间硬化模型中的t₂；通过得到主程序返回当前时刻的应变值ε，根据表达式2可迭代求得在σ₂,T₂状态下的蠕变曲线产生应变ε的对应的时间t₂′。分别用t₂和t₂′替换usercreep子程序中输出变量中的t₂，即可完成符合相对时间硬化模型和应***化模型的usercreep子程序。

3.优点功效

所发明的基于归一化参数的蠕变模型能够模拟完整的3个阶段的蠕变变形，解决当前常用的蠕变模型描述蠕变曲线的不足。通过编写usercreep子程序，能够将该模型用于实际结构的蠕变变形分析，同时针对变载情况，在子程序中实现3种变载模型，能用于计算应力松弛行为。

四.附图说明：

图1：本发明所述方法流程图；

图2：典型的蠕变曲线；

图3：归一化参数蠕变模型各部分的意义；

图4：蠕变变形表达式中各系数的意义；

图5：直接时效GH4169G的蠕变曲线(680℃下不同应力)；

图6：直接时效GH4169G的蠕变曲线(650MPa下不同温度)；

图7：用表达式1对试验曲线的拟合结果(680℃下不同应力)；

图8：用表达式1对试验曲线的拟合结果(650MPa下不同温度)；

图9：用表达式2对试验曲线的拟合结果(680℃下不同应力)；

图10：用表达式2对试验曲线的拟合结果(650MPa下不同温度)；

图11：适用于变载情况的时间硬化模型；

图12：适用于变载情况的相对时间硬化模型；

图13：适用于变载情况的应***化模型；

图中符号说明如下：

η₁，η₂，η₃为模型中参数，分别表示蠕变3个阶段的蠕变量；

ζ为无量纲时间，ζ＝t/t_c，t_c为给定温度和应力下的持久寿命，ζ∈[0,1]。

五.具体实施方案

本发明一种描述蠕变变形的方法，它是一种能完整描述3个阶段蠕变变形的方法，所述的“3个阶段”是指：第1阶段由于变形引起加工硬化导致蠕变速率随时间不断降低，称为初始蠕变阶段；第2阶段蠕变速率不随时间变化，即稳态蠕变阶段；第3阶段蠕变速率随时间增大直至断裂，称为加速蠕变阶段，其具体步骤如下：

步骤一：将试验得到的不同温度和应力下的蠕变曲线归一化；

用试验得到的不同温度和应力下的蠕变曲线的时间坐标除以该温度和应力下的持久断裂时间，则所有曲线的横坐标均为归一化时间坐标ζ＝t/t_c，ζ∈[0,1]。

步骤二：对曲线的第2阶段进行拟合，得到η₁，η₂，η₃；

对每条曲线第2阶段进行拟合，若曲线第2阶段较为明显，则得到直线方程y＝kx+b，否则，根据最小蠕变率点得到直线方程；得到的一组k值则为不同应力和温度下的一组η₂值，即归一化坐标下的稳态蠕变应变率或最小蠕变率；拟合直线中的b值即η₁，继而可得到η₃＝ε_r-η₁-η₂，ε_r为断裂时的蠕变应变。

步骤三：拟合得到η₄，η₅；

然后由已得到的参数η₁，η₂，η₃，根据表达式或对每条曲线进行拟合(可采用matlab等软件进行最小二乘拟合)，得到η₄，η₅。

步骤四：将η_i(i＝1,2,3,4,5)表示为温度和应力的函数；

将η_i(i＝1,2,3,4,5)表示为无量纲应力σ/σ_0.2与无量纲温度T/T_m的函数。

步骤五：将该方法与有限元软件结合，编写usercreep子程序；

利用通用有限元软件ANSYS中usercreep子程序，将步骤一～步骤四得到的表达式写入子程序中，以达到计算实际结构蠕变变形的目的；对所编写的子程序进行编译连接后，在ANSYS主程序中调用子程序即可用所发明的方法对实际结构进行蠕变变形计算。

步骤六：选择适当的适用于变载情况下的蠕变模型；

在usercreep子程序中实现3种不同的适用于变载情况下蠕变模型：时间硬化模型，相对时间硬化模型，应***化模型。时间硬化模型为当载荷在t时刻变化(由σ₁,T₁变到σ₂,T₂)，t时刻后的蠕变曲线由σ₂,T₂状态下t时刻后的蠕变应变曲线上下平移得到；相对时间硬化模型为当载荷在t时刻变化(由σ₁,T₁变到σ₂,T₂)，t时刻后的蠕变曲线由σ₂,T₂状态下t₂＝t×t_c,1/t_c,2时刻后的蠕变曲线平移得到，t_c,1和t_c,2分别表示σ₁,T₁和σ₂,T₂状态下的持久寿命；应***化模型为载荷在t时刻变化(由σ₁,T₁变到σ₂,T₂)，t时刻后的蠕变曲线由σ₂,T₂状态下产生与前一状态相同的蠕变应变所对应的时间之后曲线左右平移得到。

步骤七：建立实际结构的有限元模型，进行蠕变变形计算和应力松弛行为分析；

对实际结构在CAE前处理软件中建立有限元模型；由于应力松弛行为属于变载条件下的蠕变行为，能通过所编写的usercreep子程序进行计算分析；在ANSYS主程序中调用所编写的usercreep子程序，同时输入所求得模型参数值，对实际结构进行有限元数值模拟，计算蠕变变形和应力松弛行为。

其中，在步骤二中所述的“对每条曲线第2阶段进行拟合，若曲线第2阶段较为明显，则得到直线方程y＝kx+b，否则，根据最小蠕变率点得到直线方程”，其方法如下：利用第二阶段的试验数据点，可采用excel或matlab等软件对该数据点进行一次多项式拟合；若蠕变曲线的第2阶段并不明显，可采用多次多项式对蠕变曲线进行拟合，然后对多次多项式求导，导数值最小的点即为最小蠕变率点，导数值即η₂。

其中，在步骤四中所述的“将η_i(i＝1,2,3,4,5)表示为温度和应力的函数”，具体做法如下：选取特定函数将无量纲应力σ/σ_0.2与无量纲温度T/T_m表示为η_i＝f(σ/σ_0.2,T/T_m)，例如：

${\mathrm{\η}}_i=a_i+b_i\frac{T}{T_m}+c_i\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}}+d_i\frac{T}{T_m}\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}}$ 或 $\ln {\mathrm{\η}}_i=a_i+b_i\frac{T}{T_m}+c_i\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}}+d_i\frac{T}{T_m}\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}},$

其中，a_i，b_i，c_i，d_i(i＝1,2,3,4,5)为材料相关系数，可由进行最小二乘拟合(可采用matlab等软件)得到；通过一组a_i，b_i，c_i，d_i的系数值，则可求得任意温度和应力下的η_i，继而得到该温度和应力下的蠕变曲线。

其中，在步骤五中所述的“usercreep子程序”，其编写所需的输出变量为：蠕变应变增量delcr、蠕变应变增量对等效应力的导数dcrda(1)，蠕变应变增量对蠕变应变的导数dcrda(2)；表达式2(即)的各输出变量为：

蠕变应变增量delcr：

$\mathrm{delcr}={\stackrel{\·}{\mathrm{\ϵ}}}_c\mathrm{\Δt}=\frac{1}{t_c}({\mathrm{\η}}_1{\mathrm{\η}}_4{e}^{{-\mathrm{\η}}_4\mathrm{\ζ}}+{\mathrm{\η}}_2+{\mathrm{\η}}_3{\mathrm{\η}}_5{\mathrm{\ζ}}^{{\mathrm{\η}}_5-1})\mathrm{\Δt}$

蠕变应变增量对等效应力的导数dcrda(1)：

$\begin{array}{c}\mathrm{dcrda}\left(1\right)=\frac{\∂\left({\stackrel{\·}{\mathrm{\ϵ}}}_c\mathrm{\Δt}\right)}{\∂\mathrm{\σ}}\\ =\frac{1}{t_c}\frac{1}{{\mathrm{\σ}}_{0.2}}[{\mathrm{\η}}_1{\mathrm{\η}}_4{e}^{{-\mathrm{\η}}_4\mathrm{\ζ}}(q_1+q_4-q_4{\mathrm{\η}}_4\mathrm{\ζ})+q_2{\mathrm{\η}}_2+{\mathrm{\η}}_3{\mathrm{\η}}_5{\mathrm{\ζ}}^{{\mathrm{\η}}_5-1}(q_3+q_5+q_5{\mathrm{\η}}_5\ln \mathrm{\ζ})]\mathrm{\Δt}\end{array}$

式中， $q_i=c_i+d_i\frac{T}{T_m}(i=\mathrm{1,2,3,4,5});$

蠕变应变增量对蠕变应变的导数dcrda(2)：

$\mathrm{dcrda}\left(2\right)=\frac{\∂\left({\stackrel{\·}{\mathrm{\ϵ}}}_c\mathrm{\Δt}\right)}{\∂{\mathrm{\ϵ}}_c}=\frac{d{\stackrel{\·}{\mathrm{\ϵ}}}_c}{\mathrm{dt}}\left(\frac{1}{\frac{d{\mathrm{\ϵ}}_c}{\mathrm{dt}}}\right)\mathrm{\Δt}=\frac{{-\mathrm{\η}}_1{\mathrm{\η}}_4^2{e}^{{-\mathrm{\η}}_4\mathrm{\ζ}}+{\mathrm{\η}}_3{\mathrm{\η}}_5({\mathrm{\η}}_5-1){\mathrm{\ζ}}^{{\mathrm{\η}}_5-2}}{t_c({\mathrm{\η}}_1{\mathrm{\η}}_4{e}^{{-\mathrm{\η}}_4\mathrm{\ζ}}+{\mathrm{\η}}_2+{\mathrm{\η}}_3{\mathrm{\η}}_5{\mathrm{\ζ}}^{{\mathrm{\η}}_5-1})}\mathrm{\Δt}.$

其中，在步骤六中所述的“变载情况下的蠕变模型”，在子程序中实现方法为：

在子程序中，根据t₂＝t×t_c,1/t_c,2，可计算得到相对时间硬化模型中的t₂；通过得到主程序返回当前时刻的应变值ε，根据表达式2(即)，可迭代求得在σ₂,T₂状态下的蠕变曲线产生应变ε的对应的时间t₂′；分别用t₂和t₂′替换usercreep子程序中输出变量中的t₂，即可完成符合相对时间硬化模型和应***化模型的usercreep子程序。

下面将结合附图和实例对本发明作进一步的详细说明。图1为本发明所述方法的流程图。图2为典型的蠕变蠕变，蠕变3个阶段较为明显。第1阶段由于变形引起加工硬化导致蠕变速率随时间不断降低，称为初始蠕变阶段；第2阶段为直线，这是由于加工硬化与回复软化过程达到动态平衡而导致蠕变速率保持不变，即为稳态蠕变阶段；第3阶段蠕变速率随时间增大直至断裂，称为加速蠕变阶段。图3为所采用的蠕变应变表达式的各组成部分的曲线，说明其具有描述蠕变3个阶段的能力。

图4说明蠕变应变表达式各参数的物理意义，η₁，η₂，η₃，分别为蠕变3个阶段的蠕变量，η₄，η₅分别控制蠕变第1阶段和第3阶段的变化快慢，且η₅>1。步骤二拟合的直线即图中的斜线，根据各参数物理意义可得η₁，η₂，η₃的值，从而拟合得到η₄，η₅。图5和图6为直接时效GH4169G的蠕变试验曲线，图7、图8和图9、图10分别为根据技术方案步骤一到步骤四，用表达式1和表达式2描述蠕变曲线得到的结果，可看出模型曲线和试验曲线的各部分均较为接近，说明该方法能够较好地完整描述3个阶段的蠕变曲线。

图11、图12和图13分别为变载下的时间硬化模型、相对时间硬化模型和应***化模型，分别由σ₂,T₂状态下的曲线通过平移得到。

Claims (5)

1.一种描述蠕变变形的方法，其特征在于：它是一种能完整描述3个阶段蠕变变形的方法，所述的“3个阶段”是指：第1阶段由于变形引起加工硬化导致蠕变速率随时间不断降低，称为初始蠕变阶段；第2阶段蠕变速率不随时间变化，即稳态蠕变阶段；第3阶段蠕变速率随时间增大直至断裂，称为加速蠕变阶段，其具体步骤如下：

步骤一：将试验得到的不同温度和应力下的蠕变曲线归一化；

用试验得到的不同温度和应力下的蠕变曲线的时间坐标除以该温度和应力下的持久断裂时间，则所有曲线的横坐标均为归一化时间坐标ζ＝t/t_c，ζ∈[0,1]；

步骤二：对曲线的第2阶段进行拟合，得到η₁，η₂，η₃；

对每条曲线第2阶段进行拟合，若曲线第2阶段较为明显，则得到直线方程y＝kx+b，否则，根据最小蠕变率点得到直线方程；得到的一组k值则为不同应力和温度下的一组η₂值，即归一化坐标下的稳态蠕变应变率及最小蠕变率；拟合直线中的b值即η₁，继而得到η₃＝ε_r-η₁-η₂，ε_r为断裂时的蠕变应变；

步骤三：拟合得到η₄，η₅；

然后由已得到的参数η₁，η₂，η₃，根据表达式及对每条曲线进行拟合得到η₄，η₅；

步骤四：将η_i该i＝1,2,3,4,5，表示为温度和应力的函数；

将η_i该i＝1,2,3,4,5表示为无量纲应力σ/σ_0.2与无量纲温度T/T_m的函数；

步骤五：将该方法与有限元软件结合，编写usercreep子程序；

利用通用有限元软件ANSYS中usercreep子程序，将步骤一～步骤四得到的表达式写入子程序中，以达到计算实际结构蠕变变形的目的；对所编写的子程序进行编译连接后，在ANSYS主程序中调用子程序即能用所发明的方法对实际结构进行蠕变变形计算；

步骤六：选择预定的适用于变载情况下的蠕变模型；

在usercreep子程序中实现3种不同的适用于变载情况下蠕变模型：时间硬化模型，相对时间硬化模型，应***化模型；时间硬化模型为当载荷在t时刻变化(由σ₁,T₁变到σ₂,T₂)，t时刻后的蠕变曲线由σ₂,T₂状态下t时刻后的蠕变应变曲线上下平移得到；相对时间硬化模型为当载荷在t时刻变化，由σ₁,T₁变到σ₂,T₂，t时刻后的蠕变曲线由σ₂,T₂状态下t₂＝t×t_c,1/t_c,2时刻后的蠕变曲线平移得到，t_c,1和t_c,2分别表示σ₁,T₁和σ₂,T₂状态下的持久寿命；应***化模型为载荷在t时刻变化，由σ₁,T₁变到σ₂,T₂，t时刻后的蠕变曲线由σ₂,T₂状态下产生与前一状态相同的蠕变应变所对应的时间之后曲线左右平移得到；

步骤七：建立实际结构的有限元模型，进行蠕变变形计算和应力松弛行为分析；

对实际结构在CAE前处理软件中建立有限元模型；由于应力松弛行为属于变载条件下的蠕变行为，能通过所编写的usercreep子程序进行计算分析；在ANSYS主程序中调用所编写的usercreep子程序，同时输入所求得模型参数值，对实际结构进行有限元数值模拟，计算蠕变变形和应力松弛行为。

2.根据权利要求1所述的一种描述蠕变变形的方法，其特征在于：在步骤二中所述的“对每条曲线第2阶段进行拟合，若曲线第2阶段较为明显，则得到直线方程y＝kx+b，否则，根据最小蠕变率点得到直线方程”，其方法如下：利用第二阶段的试验数据点，采用excel及matlab软件对该数据点进行一次多项式拟合；若蠕变曲线的第2阶段并不明显，则采用多次多项式对蠕变曲线进行拟合，然后对多次多项式求导，导数值最小的点即为最小蠕变率点，导数值即η₂。

3.根据权利要求1所述的一种描述蠕变变形的方法，其特征在于：在步骤四中所述的“将η_i该i＝1,2,3,4,5表示为温度和应力的函数”，具体做法如下：选取特定函数将无量纲应力σ/σ_0.2与无量纲温度T/T_m表示为η_i＝f(σ/σ_0.2,T/T_m)，例如：

${\mathrm{\η}}_i=a_i+b_i\frac{T}{T_m}+c_i\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}}+d_i\frac{T}{T_m}\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}}$ 及 $\ln {\mathrm{\η}}_i=a_i+b_i\frac{T}{T_m}+c_i\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}}+d_i\frac{T}{T_m}\frac{\mathrm{\σ}}{{\mathrm{\σ}}_{0.2}},$

其中，a_i，b_i，c_i，d_i该i＝1,2,3,4,5为材料相关系数，由进行最小二乘拟合得到；通过一组a_i，b_i，c_i，d_i的系数值，则求得任意温度和应力下的η_i，继而得到该温度和应力下的蠕变曲线。

4.根据权利要求1所述的一种描述蠕变变形的方法，其特征在于：在步骤五中所述的“usercreep子程序”，其编写所需的输出变量为：蠕变应变增量delcr、蠕变应变增量对等效应力的导数dcrda(1)，蠕变应变增量对蠕变应变的导数dcrda(2)；表达式2即 ${\mathrm{\ϵ}}_c={\mathrm{\η}}_1(1-e^{-{\mathrm{\η}}_4\mathrm{\ζ}})+{\mathrm{\η}}_2\mathrm{\ζ}+{\mathrm{\η}}_3{\mathrm{\ζ}}^{{\mathrm{\η}}_5}$ 的各输出变量为：

蠕变应变增量delcr：

$\mathrm{delcr}={\stackrel{\·}{\mathrm{\ϵ}}}_c\mathrm{\Δt}=\frac{1}{t_c}({\mathrm{\η}}_1{\mathrm{\η}}_4{e}^{-{\mathrm{\η}}_4\mathrm{\ζ}}+{\mathrm{\η}}_2+{\mathrm{\η}}_3{\mathrm{\η}}_5{\mathrm{\ζ}}^{{\mathrm{\η}}_5-1})\mathrm{\Δt}$

蠕变应变增量对等效应力的导数dcrda(1)：

$\begin{array}{c}\mathrm{dcrda}=\frac{\∂\left({\stackrel{\·}{\mathrm{\ϵ}}}_c\mathrm{\Δt}\right)}{\∂\mathrm{\σ}}\\ =\frac{1}{t_c}\frac{1}{{\mathrm{\σ}}_{0.2}}[{\mathrm{\η}}_1{\mathrm{\η}}_4{e}^{-{\mathrm{\η}}_4\mathrm{\ζ}}(q_1+q_4-q_4{\mathrm{\η}}_4\mathrm{\ζ})+q_2{\mathrm{\η}}_2+{\mathrm{\η}}_3{\mathrm{\η}}_5{\mathrm{\ζ}}^{{\mathrm{\η}}_5-1}(q_3+q_5+q_5{\mathrm{\η}}_5\ln \mathrm{\ζ})]\mathrm{\Δt}\end{array}$

式中， $q_i=c_i+d_i\frac{T}{T_m}(i=\mathrm{1,2},3,4,5);$

蠕变应变增量对蠕变应变的导数dcrda(2)：

$\mathrm{dcrda}\left(2\right)=\frac{\∂\left({\stackrel{\·}{\mathrm{\ϵ}}}_c\mathrm{\Δt}\right)}{{\∂\mathrm{\ϵ}}_c}=\frac{d{\stackrel{\·}{\mathrm{\ϵ}}}_c}{\mathrm{dt}}\left(\frac{1}{\frac{{\mathrm{d\ϵ}}_c}{\mathrm{dt}}}\right)\mathrm{\Δt}=\frac{{-\mathrm{\η}}_1{\mathrm{\η}}_4^2{e}^{{-\mathrm{\η}}_4\mathrm{\ζ}}+{\mathrm{\η}}_3{\mathrm{\η}}_5({\mathrm{\η}}_5-1){\mathrm{\ζ}}^{{\mathrm{\η}}_5-2}}{t_c({\mathrm{\η}}_1{\mathrm{\η}}_4{e}^{-{\mathrm{\η}}_4\mathrm{\ζ}}+{\mathrm{\η}}_2+{\mathrm{\η}}_3{\mathrm{\η}}_5{\mathrm{\ζ}}^{{\mathrm{\η}}_5-1})}\mathrm{\Δt}.$

5.根据权利要求1所述的一种描述蠕变变形的方法，其特征在于：在步骤六中所述的“变载情况下的蠕变模型”，在子程序中实现方法为：

在子程序中，根据t₂＝t×t_c,1/t_c,2，计算得到相对时间硬化模型中的t₂；通过得到主程序返回当前时刻的应变值ε，根据表达式2即迭代求得在σ₂,T₂状态下的蠕变曲线产生应变ε的对应的时间t₂′；分别用t₂和t₂′替换usercreep子程序中输出变量中的t₂，即完成符合相对时间硬化模型和应***化模型的usercreep子程序。