一种计算工程结构的应力应变响应的方法和系统

1.本发明涉及工程结构分析技术领域，特别是涉及一种计算工程结构的应力应变响应的方法和系统。


背景技术：

2.本构理论在结构强度设计、变形分析、应力应变计算等方面不可或缺，在机械、土木、航空航天、材料、核电等学科获得了广泛的应用。目前本构模型在工程结构局部强度分析中的应用并不成熟，其应用于工程结构的应力应变计算及变形分析还面临着些许问题。
3.本构模型的积分方案分为两种：显式积分方案与隐式积分方案。显式积分方案的优点是：无需迭代，因此计算速度快、收敛性强；然而无法保证计算精度，需划分较小时间步长。隐式积分方案的优点是：通过数值迭代求解，因此能够保证计算精度，可划分较大的时间步长；然而数值迭代带来了较大的计算代价，计算速度慢且难以保证收敛性。即无论选用何种积分方案，均难以兼顾计算速度、收敛性与误差，这一问题制约着先进本构在工程结构局部应力应变计算及整体变形分析中的应用。


技术实现要素：

4.本发明的目的是提供一种计算工程结构的应力应变响应的方法和系统，兼顾了计算速度、收敛性与误差。
5.为实现上述目的，本发明提供了如下方案：
6.一种计算工程结构的应力应变响应的方法，所述应力应变响应包括所述工程结构的应力张量和应变张量；所述方法包括以下步骤：
7.确定所述工程结构在当前温度下的材料参数。
8.获取所述工程结构各状态变量的值以及各状态变量的变化率；所述状态变量用于表征所述工程结构的形变状态或热力学状态；所述状态变量的变化率用于表征所述状态变量随时间的变化程度。
9.针对任一状态变量，根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值。
10.根据所述工程结构在当前温度下的材料参数和各状态变量的值，计算得到所述工程结构的应力张量和应变张量。
11.根据各状态变量值的变化量和所述工程结构在当前温度下的材料参数，确定计算误差。
12.判断所述计算误差是否满足预设精度；若否，则调整所述预设时间步长，并跳转到步骤“针对任一状态变量，根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值”；若是，则保持所述预设时间步长不变。
13.根据屈服函数确定所述工程结构是否发生屈服；若是，则根据所述工程结构在当前温度下的材料参数和各状态变量的值，更新各状态变量的变化率；否则，保持各状态变量
的变化率不变。
14.存储各状态变量的值以及各状态变量的变化率。
15.跳转到步骤“确定所述工程结构在当前温度下的材料参数”，进行下一预设时间步长的应力应变计算，直至达到预设的计算时间。
16.可选地，所述判断所述计算误差是否满足预设精度，具体包括：
17.判断所述计算误差是否小于预设误差上限。
18.若是，则确定所述计算误差满足预设精度。
19.若否，则确定所述计算误差不满足预设精度。
20.可选地，所述判断所述计算误差是否满足预设精度，具体包括：
21.判断所述计算误差是否大于预设误差上限。
22.若是，则确定所述计算误差不满足预设精度，调低所述预设时间步长，并跳转到步骤“针对任一状态变量，根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值”。
23.否则，判断所述计算误差是否大于预设误差下限。
24.若是，则保持所述预设时间步长不变。
25.否则，调高所述预设时间步长。
26.可选地，所述状态变量包括：背应力张量、各向同性硬化标量、损伤张量和非弹性应变张量。
27.可选地，所述根据屈服函数确定所述工程结构是否发生屈服，具体包括：
28.根据应力张量偏量和背应力张量偏量，计算得到等效应力。
29.根据所述等效应力、所述各向同性硬化标量和初始屈服强度，确定所述工程结构的屈服函数值。
30.若所述屈服函数值大于等于0，则确定所述工程结构发生屈服，否则确定所述工程结构未发生屈服。
31.可选地，所述根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值，具体包括：
32.根据所述状态变量的变化率和所述预设时间步长，确定所述状态变量的变化量。
33.根据所述状态变量的值和所述状态变量的变化量，更新所述状态变量的值。
34.可选地，所述根据所述工程结构在当前温度下的材料参数和各状态变量的值，更新各状态变量的变化率，具体包括：
35.针对任一状态变量，将各状态变量的值和所述工程结构在当前温度下的材料参数，带入所述状态变量的变化率更新函数中，更新所述状态变量的变化率。
36.可选地，所述根据所述工程结构在当前温度下的材料参数和各状态变量的值，计算得到所述工程结构的应力张量和应变张量，具体包括：
37.根据所述工程结构在当前温度下的材料参数和各状态变量的值，计算得到所述工程结构的应变张量。
38.将所述工程结构在当前温度下的材料参数和所述应变张量带入应力应变关系式中，计算得到所述工程结构的应力张量。
39.可选地，所述确定所述工程结构在当前温度下的材料参数，具体包括：
40.确定若干个实采温度点；所述当前温度位于实采温度点的最小值和实采温度点的最大值之间。
41.在各实采温度点下，采集所述工程结构的材料参数。
42.根据各实采温度点下采集的所述工程结构的材料参数，确定所述工程结构在当前温度下的材料参数。
43.对应于前述的计算工程结构的应力应变响应的方法，本发明还提供了一种计算工程结构的应力应变响应的系统，所述计算工程结构的应力应变响应的系统在被计算机运行时，执行如前文所述的计算工程结构的应力应变响应的方法。
44.根据本发明提供的具体实施例，本发明公开了以下技术效果：
45.本发明提供的一种计算工程结构的应力应变响应的方法和系统，方法包括：确定工程结构在当前温度下的材料参数；获取工程结构各状态变量的值以及各状态变量的变化率；针对任一状态变量，根据状态变量的值、状态变量的变化率和预设时间步长，更新状态变量的值；根据工程结构在当前温度下的材料参数和各状态变量的值，计算得到工程结构的应力张量和应变张量；根据各状态变量值的变化量和工程结构在当前温度下的材料参数，确定计算误差；判断计算误差是否满足预设精度；若否，则调整预设时间步长，并重新计算更新状态变量的值，若是，则保持预设时间步长不变；根据屈服函数确定工程结构是否发生屈服；若是，则更新各状态变量的变化率，否则，保持各状态变量的变化率不变；存储各状态变量的值以及各状态变量的变化率；跳转到步骤“确定工程结构在当前温度下的材料参数”，进行下一预设时间步长的应力应变计算，直至达到预设的计算时间。本发明采用显式积分方案，收敛性强，计算速度快，程序编写简单，并在显式积分方案中结合了误差控制策略，在保留其计算速度快、收敛性强特点的同时，也能够保证一定的计算精度。
附图说明
46.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
47.图1为本发明实施例1提供的一种计算工程结构的应力应变响应的方法的流程图；
48.图2为本发明实施例1提供的方法中步骤s6一种实施方式的流程图；
49.图3为本发明实施例1提供的方法中步骤s6另一种实施方式的流程图；
50.图4为本发明实施例1提供的方法中粘塑性本构模型理论框架的示意图；
51.图5为本发明实施例1提供的方法中蠕变前后涡轮叶片等效应力计算结果的示意图；
52.图6为本发明实施例1提供的方法中蠕变100h后涡轮叶片等效非弹性应变计算结果的示意图；
53.图7为本发明实施例2提供的一种计算工程结构的应力应变响应的系统的结构示意图。
具体实施方式
54.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
55.本构理论在结构强度设计、变形分析、应力应变计算等方面不可或缺，在机械、土木、航空航天、材料、核电等学科获得了广泛的应用。
56.航空航天热端部件在高温复杂载荷下服役，热端零部件局部处于复杂的多轴应力状态，载荷随时间不断改变，服役时间可长达数千小时。局部应力应变场与工程结构整体变形的准确计算需要发展先进的本构理论，精确的应力应变计算结果与变形分析是结构设计、强度评价和寿命预测的基础。
57.对于航空航天零部件结构来说，不但要求设计人员充分减小零部件质量，还要求结构具有足够的强度储备以及高可靠性，任何结构的细节设计都需要准确把握，而传统的基于线弹性或经典弹塑性理论的结构设计和分析方法通常留有较大的安全裕度，难以适用于复杂结构，随着技术的发展其表现出明显不足，需不断发展先进的本构理论用于工程结构应力应变计算及变形分析。
58.精确预测结构变形需要先进材料本构支持，准确的变形分析可提高工程结构的工作效率并减小其故障概率，nasa相关研究显示：对于航空发动机，减小0.254mm的涡轮间隙，可降低1％耗油率，实现同等功率，涡轮前温度可降低10℃，以现阶段的情况来看，仅2002年即可节约16亿美元燃油费用。
59.可见先进的寿命管理技术需要本构模型提供准确的结构局部应力应变状态，先进本构模型有助于延长大修间隔与检修周期，延长零部件的服役寿命，降低工程结构的维护成本，相关研究显示：若自1988年开始使用先进寿命管理技术，可在2000年之前节约35亿美元左右的维护费用。
60.目前先进本构理论在工程结构局部强度分析中的应用并不成熟，其应用于工程结构的应力应变计算及变形分析还面临着些许问题，这些问题并没有得到较好解决，众多问题中“计算速度、收敛性与精度的平衡问题”及“材料参数的温度相关性问题”尤为突出，非弹性应力应变分析属于非线性分析，其计算速度、收敛性与精度通常难以兼顾，一般收敛性强的计算程序需要设置很小的增量步步长，极大地增大了计算代价，完成结构强度计算校核需要大量的时间；而设置较大的增量步步长通常会导致程序难以收敛，导致计算终止，无法完成强度分析，且计算精度难以保证。
61.本发明的目的是提供一种计算工程结构的应力应变响应的方法和系统，兼顾了计算速度、收敛性与误差。
62.为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
63.实施例1：
64.本实施例提供了一种计算工程结构的应力应变响应的方法，应力应变响应的计算具体包括应力张量的计算和应变张量的计算，如图1所示，该方法包括以下步骤：
65.s1、确定所述工程结构在当前温度下的材料参数。
66.s2、获取所述工程结构各状态变量的值以及各状态变量的变化率；所述状态变量用于表征所述工程结构的形变状态或热力学状态；所述状态变量的变化率用于表征所述状态变量随时间的变化程度。
67.s3、针对任一状态变量，根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值。本实施例中优选采用显式积分方案更新所述状态变量。
68.s4、根据所述工程结构在当前温度下的材料参数和各状态变量的值，计算得到所述工程结构的应力张量和应变张量。至此计算得到了当前温度下以及当前时间步长结束时工程结构所有点的应力张量和应变张量，为了计算下一时间步长的应力应变，需要确定下一时间步长开始时各状态变量的变化率。
69.s5、根据各状态变量值的变化量和所述工程结构在当前温度下的材料参数，确定计算误差。
70.s6、判断所述计算误差是否满足预设精度；若否，则执行步骤s7调整所述预设时间步长，并跳转到步骤s3“针对任一状态变量，根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值”；若是，则保持所述预设时间步长不变，直接执行步骤s8。
71.s8、根据屈服函数确定所述工程结构是否发生屈服；若是，则执行步骤s9根据所述工程结构在当前温度下的材料参数和各状态变量的值，更新各状态变量的变化率，执行步骤s10；否则，保持各状态变量的变化率不变，直接执行步骤s10。步骤s9具体包括：
72.针对任一状态变量，将所述工程结构在当前温度下的材料参数和各状态变量的值，带入所述状态变量的变化率更新函数中，更新所述状态变量的变化率。
73.s10、存储各状态变量的值以及各状态变量的变化率。随后跳转到步骤“确定所述工程结构在当前温度下的材料参数”，进行下一预设时间步长的应力应变计算，直至达到预设的计算时间。通过上述步骤可以实现对工程结构在持续计算时间内复杂温度场下的应力应变计算。
74.作为一种可选的实施方式，如图2所示的流程图，步骤s6判断所述计算误差是否满足预设精度，具体包括：
75.s61、判断所述计算误差是否小于预设误差上限。
76.若是，则执行s62确定所述计算误差满足预设精度。
77.若否，则执行s63确定所述计算误差不满足预设精度。
78.作为另一种可选的实施方式，如图3所示的流程图，步骤s6判断所述计算误差是否满足预设精度，还可以具体包括：
79.s61、判断所述计算误差是否大于预设误差上限。
80.若是，则执行步骤s62确定所述计算误差不满足预设精度，调低所述预设时间步长，并跳转到步骤s3“针对任一状态变量，根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值”。
81.否则，执行步骤s63判断所述计算误差是否大于预设误差下限。
82.若是，则执行步骤s64保持所述预设时间步长不变。
83.否则，执行步骤s65调高所述预设时间步长。
84.一些实施方式中，状态变量可以包括：背应力张量、各向同性硬化标量、损伤张量
和非弹性应变张量。
85.在步骤s6判断计算误差满足预设精度要求后，则需判断工程结构的材料是否发生屈服。若材料发生屈服，则增量步结束时状态变化率发生改变，需要对状态变量的变化率进行计算更新；若材料未发生屈服，则无需对其进行计算更新。步骤s8根据屈服函数确定所述工程结构是否发生屈服，具体包括：
86.s81、根据应力张量偏量和背应力张量偏量，计算得到等效应力。
87.s82、根据所述等效应力、所述各向同性硬化标量和初始屈服强度，确定所述工程结构的屈服函数值。
88.s83、若所述屈服函数值大于等于0，则确定所述工程结构发生屈服，否则确定所述工程结构未发生屈服。
89.步骤s3根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值，具体包括：
90.s31、根据所述状态变量的变化率和所述预设时间步长，确定所述状态变量的变化量。
91.s32、根据所述状态变量的值和所述状态变量的变化量，更新所述状态变量的值。
92.步骤s4根据所述工程结构在当前温度下的材料参数和各状态变量的值，计算得到所述工程结构的应力张量和应变张量，具体包括：
93.s41、根据所述工程结构在当前温度下的材料参数和各状态变量的值，计算得到所述工程结构的应变张量。
94.s42、将所述工程结构在当前温度下的材料参数和所述应变张量带入应力应变关系式中，计算得到所述工程结构的应力张量。
95.步骤s1确定所述工程结构在当前温度下的材料参数，具体包括：
96.s11、确定若干个实采温度点；所述当前温度位于实采温度点的最小值和实采温度点的最大值之间。
97.s12、在各实采温度点下，采集所述工程结构的材料参数。
98.s13、根据各实采温度点下采集的所述工程结构的材料参数，确定所述工程结构在当前温度下的材料参数。
99.下面以一个具体的例子，来对本实施例提供的计算工程结构的应力应变响应的方法进行说明，典型可供选择的本构模型有弹塑性本构模型、粘塑性本构模型与晶体塑性本构模型，其中弹塑性本构模型较为简单，无法描述高温结构材料时间相关的蠕变过程，晶体塑性模型过于复杂，不够通俗易懂，而粘塑性本构模型能够描述时间相关的蠕变过程，且难度适中，可较为精确地模拟热端部件局部的应力应变响应。且粘塑性本构模型可同时计算时间相关的蠕变变形与时间无关的塑性变形，在工程界获得了广泛的应用。
100.本例中以航空发动机涡轮叶片为具体实例，完整介绍粘塑性本构模型在实际工程结构中的数值实现方法，本例中计算工程结构的应力应变响应的方法主要包括以下步骤：
101.1.1材料参数对温度进行插值。
102.显式积分方案兼顾了程序的运行速度与计算精度，若要计算涡轮叶片在复杂温度场下的力学响应，仅仅依靠几组特定温度下的材料参数是远远不够的，对于“怎样利用有限的材料参数计算复杂温度场下工程结构的力学响应”问题，目前工程界并没有成熟统一的
方法。因此，在本例中还需要将材料参数对温度进行插值，将材料参数对温度进行插值的目的是获得某一温域内任意温度下的材料参数。
103.读取计算点的温度t，确定计算点温度所在的温度区间，根据温度区间两端点进行插值。若t∈[t1,t2)，t1温度下材料参数可记为t2温度下材料参数可记为t温度下的材料参数可表达为：
[0104][0105]
由此，得到了t温度下的材料参数，若计算点的温度超过最高温度，则赋予该点最高温度的材料参数，若计算点的温度低于最低温度，则赋予该点最低温度的材料参数，材料参数只进行内插，不进行外推，之后将材料参数作为粘塑性本构模型的输入，进而计算结构材料的力学响应。
[0106]
1.2读取状态变量初值及其变化率。
[0107]
将材料参数作为粘塑性本构模型的输入，将前述得到的材料参数输入粘塑性本构模型，用于计算工程结构材料的应力应变响应。应力、应变及状态变量为粘塑性本构模型的输出变量，粘塑性本构模型在计算过程中能够输出用户关心的部分状态变量，这些状态变量可用于可视化展示工程结构的非弹性应变、损伤等，因此在程序中需要对其进行计算更新，计算状态变量值之前需要进行读取其初值与变化率；现有粘塑性本构模型的基本理论框架如图4所示：
[0108]
应力应变关系表达式，定义了应力张量、应变张量之间的关系，其中σ
ij
表示应力张量，δ
ij
为kronecker张量，下标ij表示二维矩阵中的行与列(i＝1～3,j＝1～3)，下标kk为哑指标(k＝1～3)，这里采用了“指标表示法”，遵循爱因斯坦求和约定，λ为lame常数，其表达式可写为λ＝e/(1+μ)/(1-2μ)，e为弹性模量，μ为泊松比，g为剪切模量，ε
ij
表示应变张量，表示非弹性应变张量。
[0109]
等效应力表达式，将应力张量等效为一个标量，其中，σ
eq
表示等效应力，σ
ij
'表示应力张量的偏量，x表示背应力，x
ij
'表示背应力张量的偏量。
[0110]
等效非弹性应变率表达式，将非弹性应变率张量等效为一个标量，其中，表示等效非弹性应变的变化率，表示非弹性应变张量变化率。
[0111]
流动方程，定义了非弹性应变张量变化率的表达式，流动方程中包含屈服函数、背应力、各向同性硬化标量及损伤，其中，f为屈服函数，k为材料参数(拖曳应力)，n为流动方程指数，表示考虑损伤后的应力张量偏量，表示考虑损伤后的等效应力，d表示损伤。
[0112]
屈服函数，定义了屈服函数表达式，其中，r表示各向同性硬化标量，k0表示初始屈服强度。
[0113]
背应力方程，定义了背应力变化率表达式，其中，表示背应力张量的变化率，c为线性硬化项参数，a为动态硬化项参数，β为热回复项参数，r为热回复项指数，j2(x)表示背应力的第二不变量。
[0114]
各向同性硬化标量方程，定义了各向同性硬化标量变化率的表达式，其中，b为各
向同性标量速率控制参数，w为各向同性标量饱和值。
[0115]
损伤方程，分别定义了损伤变化率表达式，其中，表示考虑了损伤之后的有效应力的第二不变量，kd为损伤拖曳应力参数，nd为损伤应力指数，kd为损伤速率控制参数。
[0116]
常规粘塑性理论共设置3个内变量：背应力、各向同性硬化标量、损伤，本例中，本构计算程序共设置4个状态变量，分别为：背应力、各向同性硬化标量、损伤、非弹性应变。其中背应力主要用于描述叶片材料的运动硬化力学现象；各向同性硬化标量主要用于描述叶片材料的各向同性硬化力学现象；损伤主要用于描述叶片材料蠕变第三阶段蠕变加速现象；非弹性应变主要用于描述叶片材料的不可恢复变形。4个状态变量主要用于记录叶片材料的背应力、各向同性硬化标量、损伤和非弹性应变的大小。
[0117]
在本步骤中，读取上一增量步结束时(此增量步开始时)的背应力、各向同性硬化标量、损伤及非弹性应变状态变量及其变化率。若此增量步为初始增量步，则此处状态变量的值置为0。
[0118]
1.3更新状态变量值。
[0119]
将读取到的状态变量及其变化率带入粘塑性本构模型，可开展计算完成状态变量的计算更新。根据此增量步开始时的状态变量及其变化率，利用显式积分方案，此处采用voigt标记法，计算此增量步结束时的状态变量，共有4组状态变量，已知其在t时刻的状态变量取值及其变化率，则t+δt时刻的状态变量取值可表达为：
[0120]
背应力：其中x表示背应力，i表示背应力分量。
[0121]
各向同性硬化标量：其中r表示各向同性硬化标量。
[0122]
损伤：其中d表示损伤。
[0123]
非弹性应变：其中ε
in
表示流动方程中的非弹性应变，i表示非弹性应变分量。
[0124]
1.4组装刚度矩阵并更新应力应变。
[0125]
完成状态变量计算更新后，还需要更新应力应变才能完成计算。应力应变并不属于状态变量，计算应力应变需要根据刚度矩阵不断进行迭代。由t温度下的弹性常数(材料参数e、μ、g)，根据弹性模量e、泊松比μ及剪切模量g组装刚度矩阵，用于有限元软件迭代求解。其中，e、μ、g为粘塑性本构模型的部分材料参数，各向同性材料的刚度矩阵可表达为：
[0126][0127]
这里以铸造高温合金k417g为例，k417g为制造叶片的材料，其属于各向同性材料，仅有两个独立的弹性常数，弹性模量e与泊松比μ。
[0128]
g为剪切模量，各向同性材料的剪切模量与其弹性模量之间存在转换关系，如下式所示：
[0129][0130]
按照下式更新应力应变，这里采用了“指标标记法”，下标ij表示二维矩阵中的行与列(i＝1～3,j＝1～3)，下标kk为哑指标(k＝1～3)，遵循爱因斯坦求和约定。
[0131][0132]
1.5误差控制策略判断计算误差是否满足要求。
[0133]
完成状态变量、应力应变计算更新后，需要判定计算误差情况，若误差满足要求，可以开始下一个增量步的计算，若误差不满足要求则需要重新计算。由于采用显式积分方案，故需要使用误差控制策略，保证计算精度，在每个时间增量步中，计算误差可定义为：
[0134][0135]
其中δ表示变化量，δε
in
则表示非弹性应变的变化量，主要用于表述非弹性变形，e表示弹性模量，δj2表示应力张量的第二不变量，可将其表达为：
[0136]
其中s
ij
表示应力张量的偏量。
[0137]
应力张量可表达为：
[0138][0139]
应力张量偏量可表达为：
[0140][0141]
通过人工设定误差的上下限，进而控制显式积分方案的时间增量步长。
[0142]
用户设定的误差上下限分别记为error
low
与error
up
，误差上下限建议分别按照下式取值：
[0143]
error
low
＝1.0
×
10-5
[0144]
error
up
＝1.0
×
10-4
[0145]
若误差小于下限，则保留此步计算结果，并放大下一增量步的时间步长，将其步长置为原来步长的两倍。
[0146]
若误差大于下限且小于上限，则保留此步计算结果，保持下一增量步的时间步长不变。
[0147]
若误差大于误差上限，说明误差过大，则此步的计算结果作废，减小此步的时间步长，将此增量步的时间步长置为原来步长的0.5倍，重新进行计算。
[0148]
1.6判断是否发生屈服。
[0149]
若误差满足要求，则需判断材料是否发生屈服。若材料发生屈服，则增量步结束时
状态变化率发生改变，需要进行计算更新；若材料未发生屈服，则无需计算更新。首先需要判断材料是否发生屈服，根据增量步结束时的状态变量，判断增量步结束时材料是否屈服。屈服函数f可写表达为：
[0150][0151]
其中，r为各向同性硬化标量，k0为初始屈服强度，而等效应力可以写为：
[0152][0153]
判断材料是否屈服，若屈服函数f大于等于0则表示材料发生屈服，若屈服函数f小于0，则表示材料未屈服。若材料屈服则需要更新状态变量的变化率，并进行存储传入下一增量步，若材料未发生屈服，则无需更新状态变量的变化率，节约计算成本。
[0154]
1.7计算更新状态变量变化率。
[0155]
若上一步判断材料屈服，则根据增量步结束时(t+δt)的非弹性应变、背应力、各向同性硬化标量与损伤，计算增量步结束时的非弹性应变率、背应力率、各向同性硬化标量率与损伤率，并储存状态变量。若上一步判断材料未屈服，则结束此增量步的计算，直接进入下一增量步的计算。
[0156]
1.8储存状态变量及其变化率。
[0157]
将增量步结束时的非弹性应变、背应力、各向同性硬化标量、损伤状态的变化率存储，作为下一个增量步初始时刻的非弹性应变、背应力、各向同性硬化标量、损伤状态的变化率，供下一个增量步计算调用，并进入下一个增量步计算。
[0158]
1.9将本构模型编写成计算程序，完成了本构模型的数值实现。
[0159]
本例中通过上述1.1～1.8的步骤，完成工程结构材料在任一温度下、一个增量步之间的应力应变响应，而实际应用时，往往需要进行持续时间复杂温度场下的应力应变响应计算，因此利用商业有限元软件，结合本构程序，借助有限元软件平台进行计算，即可完成对工程结构持续时间复杂温度场下的应力应变计算。
[0160]
本例中以通用商业软件abaqus为例，计算涡轮叶片的力学响应，本构模型可以参照上述方法编写为计算程序，结合有限元软件及本构计算子程序对涡轮叶片开展强度计算，其计算结果可参考图5-图6，其中图5表示蠕变前后的等效应力分布，图6表示蠕变100h后的等效非弹性应变。
[0161]
通过上述内容，本实施例中提供的计算工程结构的应力应变响应的方法具有下列优点，将材料参数对温度进行线性插值，可利用几组材料参数计算复杂温度场下工程结构的应力应变，可将先进的本构理论更科学便捷地用于计算工程结构的力学响应；同时采用显式积分方案，收敛性强，计算速度快，程序编写简单。并在显式积分方案中结合了误差控制策略，可保证本构模型在用于结构应力应变计算时，既能够具有较快的计算速度与较高的收敛性，又能够保证一定的计算误差。能够同时兼顾计算程序的计算速度、收敛性与计算误差，程序无需划分极小的增量步，利用误差控制模块，判断误差是否满足要求，若误差不满足要求，则自动减小增量步步长，重新开始计算，直至误差满足要求，误差控制模块的设置保证了计算程序的精度；显式积分方案使计算程序结构较为简单，无需迭代计算，保证了计算速度快并且程序收敛性强。
[0162]
实施例2：
[0163]
此外，本发明实施例1的方法也可以借助于图7所示的计算工程结构的应力应变响应的系统的架构来实现。如图7所示，该计算工程结构的应力应变响应的系统可以包括材料参数确定模块、数据读取模块、状态变量更新模块、应力应变计算模块、计算误差确定模块、计算精度判断模块、屈服程度判断模块、变化率更新模块和数据存储模块；一些模块还可以有用于实现其功能的子单元，例如在应力应变计算模块中还包括应力张量计算单元和应变张量计算单元。当然，图7所示的架构只是示例性的，在实现不同的功能时，根据实际需要，可以省略图7示出的系统中的一个或至少两个组件。
[0164]
技术中的程序部分可以被认为是以可执行的代码和/或相关数据的形式而存在的“产品”或“制品”，通过计算机可读的介质所参与或实现的。有形的、永久的储存介质可以包括任何计算机、处理器、或类似设备或相关的模块所用到的内存或存储器。例如，各种半导体存储器、磁带驱动器、磁盘驱动器或者类似任何能够为软件提供存储功能的设备。
[0165]
所有软件或其中的一部分有时可能会通过网络进行通信，如互联网或其他通信网络。此类通信可以将软件从一个计算机设备或处理器加载到另一个。例如：从视频目标检测设备的一个服务器或主机计算机加载至一个计算机环境的硬件平台，或其他实现系统的计算机环境，或与提供目标检测所需要的信息相关的类似功能的系统。因此，另一种能够传递软件元素的介质也可以被用作局部设备之间的物理连接，例如光波、电波、电磁波等，通过电缆、光缆或者空气等实现传播。用来载波的物理介质如电缆、无线连接或光缆等类似设备，也可以被认为是承载软件的介质。在这里的用法除非限制了有形的“储存”介质，其他表示计算机或机器“可读介质”的术语都表示在处理器执行任何指令的过程中参与的介质。
[0166]
本文中应用了具体个例，但以上描述仅是对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；本领域的技术人员应该理解，上述本发明的各模块或各步骤可以用通用的计算机装置来实现，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。本发明不限制于任何特定的硬件和软件的结合。
[0167]
同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

技术特征：
1.一种计算工程结构的应力应变响应的方法，所述应力应变响应包括应力张量和应变张量；其特征在于，所述方法包括以下步骤：确定所述工程结构在当前温度下的材料参数；获取所述工程结构各状态变量的值以及各状态变量的变化率；所述状态变量用于表征所述工程结构的形变状态或热力学状态；所述状态变量的变化率用于表征所述状态变量随时间的变化程度；针对任一状态变量，根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值；根据所述工程结构在当前温度下的材料参数和各状态变量的值，计算得到所述工程结构的应力张量和应变张量；根据各状态变量值的变化量和所述工程结构在当前温度下的材料参数，确定计算误差；判断所述计算误差是否满足预设精度；若否，则调整所述预设时间步长，并跳转到步骤“针对任一状态变量，根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值”；若是，则保持所述预设时间步长不变；根据屈服函数确定所述工程结构是否发生屈服；若是，则根据所述工程结构在当前温度下的材料参数和各状态变量的值，更新各状态变量的变化率；否则，保持各状态变量的变化率不变；存储各状态变量的值以及各状态变量的变化率；跳转到步骤“确定所述工程结构在当前温度下的材料参数”，进行下一预设时间步长的应力应变计算，直至达到预设的计算时间。2.根据权利要求1所述的计算工程结构的应力应变响应的方法，其特征在于，所述判断所述计算误差是否满足预设精度，具体包括：判断所述计算误差是否小于预设误差上限；若是，则确定所述计算误差满足预设精度；若否，则确定所述计算误差不满足预设精度。3.根据权利要求1所述的计算工程结构的应力应变响应的方法，其特征在于，所述判断所述计算误差是否满足预设精度，具体包括：判断所述计算误差是否大于预设误差上限；若是，则确定所述计算误差不满足预设精度，调低所述预设时间步长，并跳转到步骤“针对任一状态变量，根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值”；否则，判断所述计算误差是否大于预设误差下限；若是，则保持所述预设时间步长不变；否则，调高所述预设时间步长。4.根据权利要求1所述的计算工程结构的应力应变响应的方法，其特征在于，所述状态变量包括：背应力张量、各向同性硬化标量、损伤张量和非弹性应变张量。5.根据权利要求4所述的计算工程结构的应力应变响应的方法，其特征在于，所述根据屈服函数确定所述工程结构是否发生屈服，具体包括：
根据应力张量偏量和背应力张量偏量，计算得到等效应力；根据所述等效应力、所述各向同性硬化标量和初始屈服强度，确定所述工程结构的屈服函数值；若所述屈服函数值大于等于0，则确定所述工程结构发生屈服，否则确定所述工程结构未发生屈服。6.根据权利要求1所述的计算工程结构的应力应变响应的方法，其特征在于，所述根据所述状态变量的值、所述状态变量的变化率和预设时间步长，更新所述状态变量的值，具体包括：根据所述状态变量的变化率和所述预设时间步长，确定所述状态变量的变化量；根据所述状态变量的值和所述状态变量的变化量，更新所述状态变量的值。7.根据权利要求1所述的计算工程结构的应力应变响应的方法，其特征在于，所述根据所述工程结构在当前温度下的材料参数和各状态变量的值，更新各状态变量的变化率，具体包括：针对任一状态变量，将各状态变量的值和所述工程结构在当前温度下的材料参数，带入所述状态变量的变化率更新函数中，更新所述状态变量的变化率。8.根据权利要求1所述的计算工程结构的应力应变响应的方法，其特征在于，所述根据所述工程结构在当前温度下的材料参数和各状态变量的值，计算得到所述工程结构的应力张量和应变张量，具体包括：根据所述工程结构在当前温度下的材料参数和各状态变量的值，计算得到所述工程结构的应变张量；将所述工程结构在当前温度下的材料参数和所述应变张量带入应力应变关系式中，计算得到所述工程结构的应力张量。9.根据权利要求1所述的计算工程结构的应力应变响应的方法，其特征在于，所述确定所述工程结构在当前温度下的材料参数，具体包括：确定若干个实采温度点；所述当前温度位于实采温度点的最小值和实采温度点的最大值之间；在各实采温度点下，采集所述工程结构的材料参数；根据各实采温度点下采集的所述工程结构的材料参数，确定所述工程结构在当前温度下的材料参数。10.一种计算工程结构的应力应变响应的系统，其特征在于，所述计算工程结构的应力应变响应的系统在被计算机运行时，执行如权利要求1-9任一项所述的计算工程结构的应力应变响应的方法。

技术总结
本发明提供了一种计算工程结构应力应变响应的方法和系统，涉及工程结构分析技术领域，通过将材料参数对温度进行线性插值，可利用少量的几组材料参数计算复杂温度场下工程结构应力应变；同时在显式积分方案中结合了误差控制策略，可保证本构模型在用于结构应力应变计算时，能够同时兼顾计算速度、收敛性与计算误差，程序无需划分极小的增量步，利用误差控制模块判断误差是否满足要求，若不满足要求则自动减小增量步步长，重新开始计算，直至误差满足要求，误差控制模块保证了计算程序的精度；显式积分方案使计算程序结构较为简单、无需迭代计算，保证了计算速度快并且程序收敛性强，克服了非线性本构数值计算程序收敛性与计算速度难以兼顾的问题。算速度难以兼顾的问题。算速度难以兼顾的问题。


技术研发人员：石多奇 隋天校 范永升 李振磊 杨晓光
受保护的技术使用者：北京航空航天大学
技术研发日：2023.04.12
技术公布日：2023/7/22