CN111125960B

CN111125960B - 一种gissmo材料失效模型参数优化方法

Info

Publication number: CN111125960B
Application number: CN201911368031.1A
Authority: CN
Inventors: 何恩泽; 赵清江; 周佳; 史爱民; 赵岩; 梁宾; 郭怡晖
Original assignee: China Automotive Engineering Research Institute Co Ltd
Current assignee: China Academy Of Automobile Technology Co ltd; China Automotive Engineering Research Institute Co Ltd
Priority date: 2019-12-26
Filing date: 2019-12-26
Publication date: 2022-02-22
Anticipated expiration: 2039-12-26
Also published as: CN111125960A

Abstract

本发明涉及材料技术领域，具体涉及一种GISSMO材料失效模型参数优化方法，包括步骤：S1确定仿真输入的真实应力应变曲线；S2通过单轴拉伸虚拟样件仿真结果与试验结果对比确定WF的初始范围；S3不带关键字*MAT_ADD_EROSION情况下，采用基于元模型的区间缩减序列优化材料参数WF；S4增加*MAT_ADD_EROSION即GISSMO失效模型，采用与步骤S3一致的优化方法与目标函数优化GISSMO失效模型参数。本发明基于商业有限元软件LS‑DYNA中提供的GISSMO失效模型，依据材料力学性能试验数据参数反求、校准GISSMO失效模型参数；采用LS‑OPT能快速识别材料参数，使得仿真与试验的输出工程应力应变曲线能获得较好一致性，也可以为快速自动化的高精度失效材料库的建立提供参考。

Description

一种GISSMO材料失效模型参数优化方法

技术领域

本发明涉及材料技术领域，具体涉及一种GISSMO材料失效模型参数优化方法。

背景技术

高强钢在汽车中应用比例不断扩大，但强度提高带来材料延展性降低，导致在某些碰撞工况中易发生断裂。传统耐撞性仿真分析较成熟地应用于车身结构开发，但其准确性仍然有待提高。对此，文件CN106096259B公开了一种材料失效分析方法，包括以下步骤：获取失效零件的失效背景；对失效零件进行外观检查，得到失效零件的失效信息；在材料失效分析数据库中查找与所述失效背景和/或所述失效信息匹配的材料失效数据，得到所述失效零件的预估失效原因；根据预估失效原因，对所述失效零件进行断口观察，得到所述失效零件的失效形式；采用材料分析法对所述失效零件与未失效零件进行对比分析，得到对比结果；根据所述失效形式和所述对比结果，确定所述失效零件的失效原因。利用该发明可提高材料失效分析人员的效率和水平。

但是，直接使用标准单轴拉伸的断裂延伸率来评价材料失效的传统方法会使碰撞断裂预测结果不可靠，并且在指导产品开发过程中会导致过于保守的设计。对材料失效的准确预测不仅需要不同应力三轴度下断裂曲线作为输入，而且要对GISSMO损伤模型中材料参数进行识别，而模型参数需要借助试验结果进行仿真标定，传统人工迭代调整失效参数直至计算结果与试验吻合不仅工作量大，而且不一定是最优解。

发明内容

本发明提供一种GISSMO材料失效模型参数优化方法，解决了传统借助实验结果标定模型参数时计算量大且难以获得最优解的技术问题。

本发明提供的基础方案为：一种GISSMO材料失效模型参数优化方法，包括步骤：S1确定仿真输入的真实应力应变曲线；S2通过单轴拉伸虚拟样件仿真与试验结果对比确定WF的初始范围；S3不带关键字*MAT_ADD_EROSION情况下，采用基于元模型的区间缩减序列优化材料参数WF；S4增加*MAT_ADD_EROSION即GISSMO失效模型，采用与步骤S3一致的优化方法与目标函数优化GISSMO失效模型参数。

本发明的工作原理在于：基于商业有限元软件LS-DYNA中提供的GISSMO失效模型，依据材料力学性能试验数据参数反求、校准GISSMO失效模型参数，采样点结果拟合方法采用元模型优化算法，不断迭代参数使得仿真与试验力位移曲线均方差最小。本发明的优点在于：采用LS-OPT能快速识别材料参数，使得仿真与试验的输出工程应力应变曲线能获得较好一致性，也可以为快速自动化的高精度失效材料库的建立提供参考。

本发明一种GISSMO材料失效模型参数优化方法，依据材料力学性能试验数据参数反求、校准GISSMO失效模型参数，只需要若干组实验的采样点的拟合曲线即可进行优化；而传统方法需要人工迭代，调整失效参数直至计算结果与试验吻合。相比较传统人工迭代10⁸数量级的计算量而言，本发明仅需要进行若干组实验，工作量大幅度减小；除此之外，通过曲线匹配，相对于人工迭代更容易得到最优解。

进一步，步骤S3中优化包括曲线预处理、整体搜索和局部匹配三个阶段。这样对于特征点较少的曲线，根据曲率极值点将曲线划分为多条曲线段，采用局部线性搜索法进行曲线的部分匹配；有利于减少曲线匹配的搜索区间，提高了曲线的匹配速度。

进一步，曲线预处理包括步骤：对输入的曲线离散采样；计算采样点的曲率；提取曲率极大值的点为特征点；将端点提取为特征点。这样可以快速、准确地为曲线匹配创造良好的条件。

进一步，整体搜索包括步骤：提取表曲率极大值为特征点；检查特征点，若特征点少于3个则进行局部拟合；建立特征点距离矩阵；利用子矩阵进行匹配；检验曲线匹配的均方差是否满足要求，若不满足要求结束整体搜索，若满足要求则进行局部匹配。通过提取曲线的特征点从而获得曲线的整体信息，用力与保证匹配的准确性。

进一步，局部匹配包括步骤：候选区域等弧长采样；计算采样点曲率；比较曲率差的平方和；扩展与对齐。局部匹配利用曲率极大值点将曲线划分为不包含极大值点的子线段，排除了包含曲率极大值点的搜索区间，减少了计算量。

进一步，步骤S1中的仿真曲线包括紧缩前半段和延拓后半段，紧缩前半段为金属材料拉伸试验获得的数据处理成的真实应力应变曲线，延拓后半段在紧缩前半段基础上延拓而来。通过定义延拓段上边界是失稳前试验真实应力应变曲线最后两点切线延长线，延拓段下边界是使用Ludwik power law硬化方程进行数据外延，上下边界应变应值均外推到1。这样实际仿真输入的真实应力由试验获得的紧缩前半段以及紧缩后半段上下边界加权平均而来。

进一步，步骤S3使用线性多项式元模型以及D-optimal采集样本点。这样采集样本点可以防止震荡和无法收敛的问题。

进一步，步骤S3优化的方法是使得试验目标工程应力应变曲线与仿真输出的工程应力应变曲线匹配度或均方差达到预设的要求。这样的优化方法可以使得试验目标工程应力应变曲线与仿真输出的工程应力应变曲线匹配度越高，即均方差最小。

进一步，步骤S4优化的参数包括应力衰减指数FADEXP、断裂曲线的尺寸效应系数SFO和尺寸效应曲线。尺寸效应曲线，横坐标为网格尺寸大小、纵坐标为对应网格修正系数，该曲线能够人为修正因网格尺寸大小不同带来的不真实的损伤累计，从而准确地预料材料的失效。

进一步，还包括步骤S5，将优化的参数和尺寸效应曲线以Excel表格的形式输出。Excel形式的文件便于读取，也方便观察、对比结果。

附图说明

图1为本发明一种GISSMO材料失效模型参数优化方法实施例的流程图。

具体实施方式

下面通过具体实施方式进一步详细的说明。

实施例1

金属材料断裂失效机理可分为正向断裂、剪切断裂以及失稳失效。正向断裂是金属变形过程中内部孔洞的形成并联合导致的，其断口形貌较为粗糙；而剪切断裂是由于剪切带滑移引起，其断口形貌较为光滑。GISSMO失效模型对材料断裂行为提供了一个唯象方程的描述，假设在受力的截面为S区域中包含部分微缺陷，截面实际的有效区域为

在本构方程中引入了损伤参数D，则损伤系数可表示为：

有效区域截面的减少带来了截面刚度的稀释，实际的应力值修正为：

σ^*＝σ(1-D) (2)

GISSMO失效模型从材料屈服、硬化、失稳直至破坏的过程中描述各种金属材料损伤累积行为。在数值仿真软件LS-DYNA中通过定义关键字*MAT_ADD_EROSION与*MAT_024、*MAT_036、*MAT_102弹塑性模型耦合使用，分别描述材料的弹塑性行为以及断裂行为，为材料损伤描述提供最大的可能性。该模型基于以下损伤累积的增量公式：

式中ε_f为与应力三轴度有关的材料断裂时刻等效塑性应变，Δε_f是每个时间步长内等效塑性应变增量，D为材料损伤，ΔD为每个时间步长内损伤变量增量，DMGEXP为非线性损伤累积指数。

此外，GISSMO失效模型一个重要的方面在于它能对材料失稳参数进行累积，用字母F表示。通过失稳参数进一步地对真实应力曲线修正，失稳参数累积的增量公式如下：

式中，ΔF是每个时间步长内失稳参数增量，F是失稳参数，DMGEXP是非线性损伤累积指数，ε_p,loc是与应力三轴度相关的材料失稳时刻等效塑性应变，Δε_p是每个时间步长内等效塑性应变增量。

失稳参数F根据式子(4)进行累积得到。当F＝1时，GISMMO模型开始对真实应力进行修正。当等效塑性应变达到材料失稳临界塑性应变时，流动应力会因为损伤的累积而减小。应力修正值与失稳损伤之间的耦合关系公式如下：

式中

为原始等效应力，σ为修正等效应力，FADEXP为应力衰减指数，D_CRIT是当失稳参数F累积到1时刻的损伤变量。令FADEXP＝1，DCRIT＝0，可得勒梅特方程，如下式

确定仿真输入的真实应力应变曲线，该曲线由前半段和后半段两段构成。前半段来自GB/T 228.1-2010金属材料拉伸试验获得的目标材料静态拉伸曲线，后半段则是在前半段曲线基础上进行延拓的延拓段。延拓段的上边界是失稳前试验真实应力应变曲线最后两点切线延长线，下边界是使用Ludwik power law硬化方程进行数据外延得到，上、下边界应变应值均外推到1。实际仿真输入的真实应力由试验获得的紧缩前半段以及紧缩后半段上下边界加权平均而来，加权系数为WF，则标定曲线＝下边界曲线*(1-WF)+上边界*WF。

通过单轴虚拟样件仿真大概确定WF的取值范围。由于WF的取值变化区间在试验工程应力应变曲线上下波动，这样可以缩小该参数优化范围，降低运算量，也可以减少该参数优化范围。

在不带关键字*MAT_ADD_EROSION情况下，采用基于元模型的区间缩减序列优化策略，使用线性多项式元模型以及D-optimal采集样本点优化材料参数WF。采样点的结果拟合采用元模型(Metamodel-based optimization)，也称响应面，是指输出响应变量Y与一组输入变量(X1，X2…Xn)之间的函数关系。前次迭代的采样点被忽略，区间缩减序列优化设计变量区间在每个迭代步缩减，当次所有采样点都在子区间内。

优化材料参数WF的目标是使得试验目标工程应力应变曲线与仿真输出的工程应力应变曲线的匹配度或均方差达到预设的要求，尽可能达到匹配度最高、均方差最小。优化包括整体搜索和局部匹配，包括取下预处理、整体搜索和局部匹配三个阶段。首先利用曲线特征点间的距离矩阵，确定候选的匹配区间，然后通过比较曲线段的曲率进行匹配，最后根据匹配的对应点集计算交换矩阵。这样有利于减少曲线匹配的搜索区间，提高了曲线的匹配速度.

整体搜索包括：步骤一，提取表曲率极大值为特征点；步骤二，检查特征点，若特征点少于3个则进行局部拟合；步骤三，建立特征点距离矩阵；步骤四，利用子矩阵进行匹配；步骤五，检验曲线匹配的均方差是否满足要求，若不满足要求结束整体搜索，若满足要求则进行局部匹配。通过提取曲线的特征点从而获得曲线的整体信息，用力与保证匹配的准确性。

局部匹配包括，步骤一，候选区域等弧长采样；步骤二，计算采样点曲率；步骤三，比较曲率差的平方和；步骤四，扩展与对齐。相对于整体搜索，局部匹配利用曲率极大值点将曲线划分为不包含极大值点的子线段，排除了包含曲率极大值点的搜索区间，减少了计算量。

在确定真应力应变曲线失效段权重系数WF系数的基础上增加GISSMO失效模型，优化应力衰减指数FADEXP、断裂曲线的尺寸效应系数SFO、尺寸效应曲线。优化完毕后，将优化的参数和尺寸效应曲线以Excel表格的形式输出。

FADEXP会对材料工程应力应变曲线紧缩后端趋势产生影响，应力衰减指数FADEXP越大，流动应力衰减越慢，且大到一定程度FADEXP由于衰减效应越慢。

GISSMO损伤模型中单元的尺寸效应，即在紧缩后由于网格划分尺寸不同，导致累积损伤结果对表现出网格尺寸的依耐性。内在原因由于没有尺寸效应曲线修正，紧缩后同一时刻同一区域小尺寸划分单元塑性应变大于大尺寸单元划分的。因此小尺寸单元划分危险区域在数值仿真中累积损伤率先达到1而失效，宏观表现在工程应力应变曲线上小尺寸划分的样件延伸率小于大尺寸网格的。在使用场景上，由于GISSMO失效模型中参考试样尺寸所用单元尺寸为0.5mm，而在实际整车碰撞过程中，划分单元尺寸要大得多，从2.5mm到10mm不等。这就要求在需要输入一条曲线人工校准损伤结果对网格尺寸的依耐性，而保障不同尺寸划分区域的损伤值的统一性。

基于以上原因，为了人工校准网格尺寸影响，使用了5个不同划分尺寸的A50拉伸样件，平行段划分尺寸分别为0.5mm、1.0mm、2.0mm、3.2mm、6.25mm。关键字中提供了卡片LCREGD，输入二维曲线，纵坐标为系数

式中ε_0.5mm为0.5mm单元划分单轴拉伸样件第一个失效单元的等效塑性应变，ε_mesh为特定尺寸划分单轴拉伸样件第一个失效单元的等效塑性应变。

通常延性金属单轴拉伸分为3个阶段，第一个是弹性变形阶段、屈服变形、均匀塑性变形阶段、局部集中塑性变形。从屈服点开始到塑性失稳点，即在均匀塑性变形阶段，塑性失稳点的应力应变仍然可以用上述公式求出，工程应力应变根据体积不变架势只能转换到样件发生失稳之前即均塑性变形。真实应力、应变曲线只能描述到紧缩之前，紧缩发生后，标距段的变形不再均匀，紧缩发生之前塑形变形在标距段是均匀的变化即满足体积不变假设，紧缩开始后为集中塑形变形阶段，不满足体积不变假设。

为获得集中塑性变形阶段材料真实应力应变曲线，需采取数据外延方法对已有曲线进行外拓。通常有三种方法对曲线进行外拓：(1)凭工程师经验根据低应变下测得试验曲线趋势进行外延；(2)使用数值分析的方法根据前面测得曲线进行拟合得到拟合公式，通过公式进行外延；(3)将测得数据代入某一材料本构方程，求得本构方程常数，再通过本构方程来求去后段曲线。为获得集中塑性变形后端曲线上下边界，上边界是失稳前真实应力应变曲线最后两点切线延长至1，下边界是使用Ludwik power law进行数据外延，应变外推到1。

对于板料，加工硬化过程通常可以使用幂指数方程定义：即用Ludwik Power Law数学描述，

K是强化系数n是加工硬化指数。

在不带GISSMO材料模型情况下，WF越大曲线后端平均斜率越大，意味着硬化效应明显，同样的变形需要更大的屈服去抵抗变形。

实施例2

与实施例1相比，不同之处仅在于，当采用比例加载方式时，DMGEXP增大非线性累积是一个凹函数，在断裂之前DMGEXP越大损伤累积前端慢后端块。当DMGEXP分别取2.0、2.2、2.4、2.6时，工程应力应变曲线随着DMGEXP增大在断裂稍微延迟。DMGEXP增大在体现在工程应力应变曲线下掉稍微延迟，当DMGEXP＝2时损伤非线性累积与试验结果较为符合。为减少优化模型非线性程度，DMGEXP可取2。

实施例3

与实施例2不同之处仅在于：采用元模型(Metamodel-based optimization)拟合采样点，输出响应变量Y与一组输入变量(X1，X2…Xn)之间的函数关系时，不直接输曲线；而是先拟合出两条曲线Lm和LM，然后Lm和LM之间形成一个区域，在该区域之间挑出某条曲线LR，并将该曲线LR作为试验目标工程应力应变曲线。

由于测得的采样点数据会有一定的误差，实际的采样点的数值必然位于某个数值区间Xim≤Xi≤XiM，Yim≤Yi≤YiM。也即Xi∈【Xim，XiM】，Xim为实验时采集到的Xi的最小值，XiM为实验时采集到的Xi的最大值；Yi∈【Yim，YiM】，Yim为实验时采集到的Yi的最小值，YiM为实验时采集到的Yi的最大值；i＝1、2、3…n。因此，若直接将这些采样点拟合成一条曲线，并将该曲线作为试验目标工程应力应变曲线，必然会产生误差。

为降低由直接拟合造成的误差，拟合采样点时按照以下步骤进行：Sa，将实验时采集到的Xi的最小值Xim拟合成曲线Lm；Sb，将实验时采集到的Xi的最大值XiM拟合成曲线LM；Sc，在Lm和LM形成的区域之间挑出某条曲线LR，并将该曲线LR作为试验目标工程应力应变曲线。比如，挑选时采用最简单的加权法，即LR＝βLm+(1-β)×LM，β为系数，且0＜β＜1。可见，Lm＜LR＜LM，这样得出的LR曲线，位于Lm和LM之间的某个位置，将其作为试验目标工程应力应变曲线可以降低误差。

除此之外，还能确定出LR曲线位于Lm和LM曲线之间的概率P。由于在Lm和LM曲线附近必然存在某些未被拟合的采样点，就第i个采样点而言，该采样点对应处Lm和LM之间的区间【Lmi，LMi】一定在区间【Yim，YiM】之内，即【Lmi，LMi】∈【Yim，YiM】。可见，若Yim和YiM未被拟合在曲线Lm和LM上，区间【Lmi，LMi】与区间【Yim，YiM】必定不重合。Yim和YiM被拟合，存在某个概率Pi。对于概率的概率，可以通过贝塔分布估计，即Pi～Be(α1i，α2i)，α1i、α2i为参数，α1i＝1，α2i＝1。

首先，计算Yim被拟合在曲线Lm上的概率PiA。第一步，分别收集第i处Lm曲线附近的采样点数量N1i；第二步，收集第i处采样点落在Lm曲线上的数量N2i；第三步，生成对应的贝塔函数，这时α1i＝N2i+1、α2i＝N1i-N2i+1，即PiA～Be(N2i+1，N1i-N2i+1)；第四步，根据贝塔分布表格计算出概率PiA。

然后，类似的，根据以上步骤，重复第一步到第四步，计算YiM被拟合在曲线LM上的概率PiB。最后可得Yim和YiM同时被拟合的概率Pi＝PiA×PiB，Pi也即区间【Lmi，LMi】位于区间【Yim，YiM】的概率。

由于各采样点相互独立，考虑到i＝1、2、3…n，可得LR曲线位于Lm和LM曲线之间的概率P＝P1×P2×P3×…Pn，即P＝∏Pi。

以上所述的仅是本发明的实施例，方案中公知的具体结构及特性等常识在此未作过多描述，所属领域普通技术人员知晓申请日或者优先权日之前发明所属技术领域所有的普通技术知识，能够获知该领域中所有的现有技术，并且具有应用该日期之前常规实验手段的能力，所属领域普通技术人员可以在本申请给出的启示下，结合自身能力完善并实施本方案，一些典型的公知结构或者公知方法不应当成为所属领域普通技术人员实施本申请的障碍。应当指出，对于本领域的技术人员来说，在不脱离本发明结构的前提下，还可以作出若干变形和改进，这些也应该视为本发明的保护范围，这些都不会影响本发明实施的效果和专利的实用性。本申请要求的保护范围应当以其权利要求的内容为准，说明书中的具体实施方式等记载可以用于解释权利要求的内容。

Claims

1.一种GISSMO材料失效模型参数优化方法，包括步骤：

S1确定仿真输入的真实应力应变曲线；

S2通过单轴拉伸虚拟样件仿真与试验结果对比确定WF的初始范围；

S3不带关键字*MAT_ADD_EROSION情况下，采用基于元模型的区间缩减序列优化材料参数WF；

S4增加*MAT_ADD_EROSION即GISSMO失效模型，采用与步骤S3一致的优化方法与目标函数优化GISSMO失效模型参数；

其特征在于，S3中，采用元模型拟合采样点，先拟合出两条曲线Lm和LM，在Lm和LM之间形成一个区域，在该区域之间挑出曲线LR，将曲线LR作为试验目标工程应力应变曲线，具体包括：

Sa，将实验时采集到的Xi的最小值Xim拟合成曲线Lm；

Sb，将实验时采集到的Xi的最大值XiM拟合成曲线LM；

Sc，在Lm和LM形成的区域之间挑出某条曲线LR，LR＝βLm+(1-β)×LM，β为系数，且0＜β＜1；Xi为第i个输入变量的值，i=1、2、3…n，n为输入变量的总个数。

2.根据权利要求1所述的一种GISSMO材料失效模型参数优化方法，其特征在于：步骤S3中优化包括曲线预处理、整体搜索和局部匹配三个阶段。

3.根据权利要求2所述的一种GISSMO材料失效模型参数优化方法，其特征在于：曲线预处理包括步骤：对输入的曲线离散采样；计算采样点的曲率；提取曲率极大值的点为特征点；将端点提取为特征点。

4.根据权利要求3所述的一种GISSMO材料失效模型参数优化方法，其特征在于：整体搜索包括步骤：提取表曲率极大值为特征点；检查特征点，若特征点少于3个则进行局部拟合；建立特征点距离矩阵；利用子矩阵进行匹配；检验曲线匹配的均方差是否满足要求，若不满足要求结束整体搜索，若满足要求则进行局部匹配。

5.根据权利要求4所述的一种GISSMO材料失效模型参数优化方法，其特征在于：局部匹配包括步骤：候选区域等弧长采样；计算采样点曲率；比较曲率差的平方和；扩展与对齐。

6.根据权利要求5所述的一种GISSMO材料失效模型参数优化方法，其特征在于：步骤S1中的仿真曲线包括紧缩前半段和延拓后半段，紧缩前半段为金属材料拉伸试验获得的数据处理成的真实应力应变曲线，延拓后半段在紧缩前半段基础上延拓而来。

7.根据权利要求6所述的一种GISSMO材料失效模型参数优化方法，其特征在于：步骤S3使用线性多项式元模型以及D-optimal采集样本点。

8.根据权利要求7所述的一种GISSMO材料失效模型参数优化方法，其特征在于：步骤S3优化的方法是使得试验目标工程应力应变曲线与仿真输出的工程应力应变曲线匹配度或均方差达到预设的要求。

9.根据权利要求8所述的一种GISSMO材料失效模型参数优化方法，其特征在于：步骤S4优化的参数包括应力衰减指数FADEXP、断裂曲线的尺寸效应系数SFO和尺寸效应曲线。

10.根据权利要求9所述的一种GISSMO材料失效模型参数优化方法，其特征在于：还包括步骤S5，将优化的参数和尺寸效应曲线以Excel表格的形式输出。