加载历史相关的数据驱动多尺度并发模拟方法及系统

1.本发明涉及材料力学行为分析技术领域，具体地，涉及一种加载历史相关的数据驱动多尺度并发模拟方法及系统。


背景技术：

2.随着材料科学的发展，材料的微观组织越来越精细，微观组织同时也决定着材料的强化机制和宏观性能。例如对于纤维增强复合材料而言，在不同尺度上存在着树脂交联网状结构、单向纤维增强结构、纤维编织结构以及零部件几何结构等；对于铝合金材料而言，在不同尺度上存在着固溶强化、颗粒强化、形变强化、细晶强化等多种强化机制。材料的多尺度模拟方法能够在不同尺度间建立材料的组织和性能的内在联系，成为近年来的研究热点。代表性体积元方法(representat ive volume element,rve)是多尺度分析中的一种常用方法，通过建立一个代表材料局部细观特征的周期性有限元模型，来预测材料在细观尺度上的力学响应演化过程。然而，在多尺度并发计算时，需要在宏观尺度每个增量步的每个积分点处调用rve模型的均质化本构参与计算。传统的高保真rve模型计算量太大，求解成本太高，无法满足多尺度并发计算的需求，迫切需要对rve进行降阶处理，以提升其求解速度。
3.自适应聚类分析(self-consistent clustering analysis,sca)方法是一种高精高效的材料多尺度并发模拟方法，能够显著提高rve的计算效率，有效降低求解成本，同时具有良好的求解精度。在该方法中，首先通过机器学习的聚类算法对rve划分集群，同一个集群内认为具有均匀的应力应变，从而降低模型的自由度数量；之后结合宏观积分点处的均质化应变约束条件，通过求解lippmann-schwinger方程，获得微观各个集群的应力应变分布；最后根据集群应力应变分布，得到rve的均质化应力，并将其返回宏观积分点，实现材料的宏微观多尺度并发模拟。自适应聚类分析方法能够从多个尺度深入分析材料内部各组成相复杂的弹塑性、损伤等非线性力学行为。然而，在自适应聚类分析过程中，仍然需要提供rve各组成相材料的唯象本构模型。随着材料的力学行为越来越复杂，现有的唯象本构模型往往难以精确描述材料实际力学行为。另外，唯象本构模型一般具有多个模型参数需要通过实验进行标定，参数标定过程又会带来新的偏差，这些偏差会在多尺度并发模拟过程中进一步放大，影响分析结果的准确性和可靠性。
4.数据驱动的计算力学方法是一种新兴的不依赖与材料本构模型的数值求解范式，在近年来得到快速的发展和完善。该方法直接从实验所获得的材料数据集出发，结合相容性条件和平衡方程等约束条件，通过对系统自由能的最小化，获取系统平衡状态的最优解。该方法无需建立材料的唯象本构模型，可直接通过材料的离散应力应变数据点求解结构的力学响应，能够从根本上避免唯象本构模型及材料参数标定所带来的系统偏差。然而，作为一种新型计算范式，数据驱动计算力学方法体系还不及有限元法成熟，在处理材料塑性等加载历史相关的力学行为时，该方法还存在一定的局限性；另一方面，目前的数据驱动计算力学方法理论体系仅局限与单一尺度下的结构仿真分析，尚无法应用于材料的多尺度并发
模拟。因此，建立一种加载历史相关的数据驱动多尺度并发模拟方法，具有十分重要的意义。


技术实现要素：

5.针对现有技术中的缺陷，本发明提供一种加载历史相关的数据驱动多尺度并发模拟方法及系统。
6.根据本发明提供的一种加载历史相关的数据驱动多尺度并发模拟方法及系统，所述方案如下：
7.第一方面，提供了一种加载历史相关的数据驱动多尺度并发模拟方法，所述方法包括：线下阶段和线上阶段；
8.其中，线下阶段包括：
9.步骤s1.1：建立代表材料组织细观特征的rve模型，计算其弹性力学响应；
10.步骤s1.2：根据rve的弹性应变分布，计算积分点处的应变集中张量；
11.步骤s1.3：基于应变集中张量，通过k-means聚类算法对积分点划分集群；
12.步骤s1.4：通过格林函数和傅里叶变换方法，计算集群相互作用张量分量；
13.线上阶段包括：
14.步骤s2.1：输入初始参考材料lame常数和线下阶段生成的相互作用张量分量；
15.步骤s2.2：根据参考材料lame常数以及线下阶段生成的相互作用张量分量，计算各个集群之间的相互作用张量；
16.步骤s2.3：通过数据驱动的方法求解lippmann-schwinger方程；
17.步骤s2.4：根据t+1时刻的集群平衡约束点计算rve均质化lame常数；
18.步骤s2.5：根据rve均质化lame常数更新参考材料lame常数，重复步骤s2.2～步骤s2.4，直到rve均质化lame常数与参考材料lame常数保持一致，得到t+1时刻最终的集群平衡约束点和rve均质化应力；
19.步骤s2.6：将rve的均质化应力返回宏观积分点，获得材料的宏微观力学响应。
20.优选地，所述步骤s1.2包括：根据六种正交加载条件下的应变分布εm(x)，计算积分点处的应变集中张量a(x)；其中应变集中张量a(x)是联系宏观均质化应变εm和微观应变分布εm(x)的四阶张量，通过以下表达式计算：
21.εm(x)＝a(x):εm。
22.优选地，所述步骤s1.3包括：
23.基于积分点的应变集中张量a(x)，通过k-means聚类算法对rve中的积分点划分集群，即完成以下关于应变集中张量的最小二乘优化过程：
[0024][0025]
其中，s代表最后的集群划分结果，sj代表第j个集群，an代表第n个积分点的应变集中张量，代表第j个集群内所有积分点的平均应变集中张量，n为总共的集群数量，||a||代表任意二阶张量a的2范数，即a的各个元素平方和的1/2次方。
[0026]
优选地，所述步骤s1.4包括：
[0027]
结合格林函数和傅里叶变换方法，计算各个集群之间的相互作用张量分量和其表达式分别为：
[0028][0029][0030]
其中，ω为rve的体积，ci为第i个集群的体积分数，χi(x)为第i个集群的特征函数，φ1和φ2为格林函数中与参考材料本构无关的部分，在傅里叶空间中具有以下形式：
[0031][0032][0033]
其中，ξ为傅里叶空间中的坐标向量，ξi是傅里叶空间中第i个方向的坐标，|ξ|为傅里叶空间坐标距离原点的距离，δ
ij
为kronecker delta符号，满足
[0034][0035]
优选地，所述步骤s2.3包括：
[0036]
步骤s2.3.1：根据宏观积分点的应变约束条件，输入微观尺度rve的均质化应变增量；
[0037]
步骤s2.3.2：根据初始材料数据集、初始材料状态变量以及t时刻集群材料状态变量，更新t时刻集群材料数据集；
[0038]
步骤s2.3.3：根据t时刻集群材料数据集，为各个集群分配初始的材料数据点；
[0039]
步骤s2.3.4：以集群应力应变状态与数据点应力应变状态的能量偏差最小作为优化目标，以lippmann-schwinger方程和均质化应变约束作为约束条件，建立约束条件下的目标函数；
[0040]
步骤s2.3.5：利用拉格朗日乘子法求解上述约束条件下的目标函数最小化问题，通过对目标函数求极小值点，建立关于集群应力应变状态的线性方程组；
[0041]
步骤s2.3.6：求解以上线性方程组，获得满足约束条件的集群应力应变分布状态，并将其作为新的集群平衡约束点；
[0042]
步骤s2.3.7：根据集群平衡约束点，在t时刻集群材料数据集中找到能量偏差最近的数据点，并将其作为新的集群材料数据点；
[0043]
步骤s2.3.8：重复步骤s2.3.6～步骤s2.3.7，直到集群的材料数据点达到收敛，得到t+1时刻的集群材料数据点和集群平衡约束点；
[0044]
步骤s2.3.9：根据t+1时刻的集群材料数据点更新集群材料状态变量；
[0045]
步骤s2.3.10：根据t+1时刻的集群平衡约束点更新rve均质化应力。
[0046]
第二方面，提供了一种加载历史相关的数据驱动多尺度并发模拟系统，所述系统
包括：线下阶段和线上阶段；
[0047]
其中，线下阶段包括：
[0048]
模块m1.1：建立代表材料组织细观特征的rve模型，计算其弹性力学响应；
[0049]
模块m1.2：根据rve的弹性应变分布，计算积分点处的应变集中张量；
[0050]
模块m1.3：基于应变集中张量，通过k-means聚类算法对积分点划分集群；
[0051]
模块m1.4：通过格林函数和傅里叶变换方法，计算集群相互作用张量分量；
[0052]
线上阶段包括：
[0053]
模块m2.1：输入初始参考材料lame常数和线下阶段生成的相互作用张量分量；
[0054]
模块m2.2：根据参考材料lame常数以及线下阶段生成的相互作用张量分量，计算各个集群之间的相互作用张量；
[0055]
模块m2.3：通过数据驱动的方法求解lippmann-schwinger方程；
[0056]
模块m2.4：根据t+1时刻的集群平衡约束点计算rve均质化lame常数；
[0057]
模块m2.5：根据rve均质化lame常数更新参考材料lame常数，重复模块m2.2～模块m2.4，直到rve均质化lame常数与参考材料lame常数保持一致，得到t+1时刻最终的集群平衡约束点和rve均质化应力；
[0058]
模块m2.6：将rve的均质化应力返回宏观积分点，获得材料的宏微观力学响应。
[0059]
优选地，所述模块m1.2包括：根据六种正交加载条件下的应变分布εm(x)，计算积分点处的应变集中张量a(x)；其中应变集中张量a(x)是联系宏观均质化应变εm和微观应变分布εm(x)的四阶张量，通过以下表达式计算：
[0060]
εm(x)＝a(x):εm。
[0061]
优选地，所述模块m1.3包括：
[0062]
基于积分点的应变集中张量a(x)，通过k-means聚类算法对rve中的积分点划分集群，即完成以下关于应变集中张量的最小二乘优化过程：
[0063][0064]
其中，s代表最后的集群划分结果，sj代表第j个集群，an代表第n个积分点的应变集中张量，代表第j个集群内所有积分点的平均应变集中张量，n为总共的集群数量，||a||代表任意二阶张量a的2范数，即a的各个元素平方和的1/2次方。
[0065]
优选地，所述模块m1.4包括：
[0066]
结合格林函数和傅里叶变换方法，计算各个集群之间的相互作用张量分量和其表达式分别为：
[0067][0068][0069]
其中，ω为rve的体积，ci为第i个集群的体积分数，χi(x)为第i个集群的特征函数，φ1和φ2为格林函数中与参考材料本构无关的部分，在傅里叶空间中具有以下形式：
[0070][0071][0072]
其中，ξ为傅里叶空间中的坐标向量，ξi是傅里叶空间中第i个方向的坐标，|ξ|为傅里叶空间坐标距离原点的距离，δ
ij
为kronecker delta符号，满足
[0073][0074]
优选地，所述模块m2.3包括：
[0075]
模块m2.3.1：根据宏观积分点的应变约束条件，输入微观尺度rve的均质化应变增量；
[0076]
模块m2.3.2：根据初始材料数据集、初始材料状态变量以及t时刻集群材料状态变量，更新t时刻集群材料数据集；
[0077]
模块m2.3.3：根据t时刻集群材料数据集，为各个集群分配初始的材料数据点；
[0078]
模块m2.3.4：以集群应力应变状态与数据点应力应变状态的能量偏差最小作为优化目标，以lippmann-schwinger方程和均质化应变约束作为约束条件，建立约束条件下的目标函数；
[0079]
模块m2.3.5：利用拉格朗日乘子法求解上述约束条件下的目标函数最小化问题，通过对目标函数求极小值点，建立关于集群应力应变状态的线性方程组；
[0080]
模块m2.3.6：求解以上线性方程组，获得满足约束条件的集群应力应变分布状态，并将其作为新的集群平衡约束点；
[0081]
模块m2.3.7：根据集群平衡约束点，在t时刻集群材料数据集中找到能量偏差最近的数据点，并将其作为新的集群材料数据点；
[0082]
模块m2.3.8：重复模块m2.3.6～模块m2.3.7，直到集群的材料数据点达到收敛，得到t+1时刻的集群材料数据点和集群平衡约束点；
[0083]
模块m2.3.9：根据t+1时刻的集群材料数据点更新集群材料状态变量；
[0084]
模块m2.3.10：根据t+1时刻的集群平衡约束点更新rve均质化应力。
[0085]
与现有技术相比，本发明具有如下的有益效果：
[0086]
1、本发明能够显著提升rve的求解效率，有效降低求解成本，在现有的有限元平台实现材料的宏微观多尺度并发模拟；
[0087]
2、本发明可以直接通过一系列离散的数据点，预测材料的宏微观力学响应，无需任何唯象本构模型，能够完全避免多尺度分析中材料本构模型及参数标定过程所导致的系统偏差；
[0088]
3、本发明可以预测加载过程中的材料塑性变形及硬化等这一类与加载历史相关的复杂力学行为，进一步扩展了数据驱动方法的应用范围和领域。
附图说明
[0089]
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、
目的和优点将会变得更明显：
[0090]
图1为加载历史相关的数据驱动多尺度并发模拟方法线下阶段流程图；
[0091]
图2为加载历史相关的数据驱动多尺度并发模拟方法线上阶段流程图；
[0092]
图3为线上阶段中通过数据驱动方法求解lippmann-schwinger方程步骤流程图；
[0093]
图4为多尺度并发模拟过程中的宏观有限元模型；
[0094]
图5为多尺度并发模拟过程中的微观rve模型；
[0095]
图6为循环加载过程中加载点的位移-时间曲线；
[0096]
图7为循环加载最终时刻宏观有限元模型；
[0097]
图8为循环加载最终时刻微观rve模型中沿x方向的应力分布；
[0098]
图9为循环加载过程中基体的应力-应变曲线；
[0099]
图10为循环加载过程中纤维的应力-应变曲线。
具体实施方式
[0100]
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
[0101]
本发明实施例提供了一种加载历史相关的数据驱动多尺度并发模拟方法，该方法包括：线下阶段和线上阶段，具体如下：
[0102]
参照图1所示，线下阶段包括：
[0103]
步骤s1.1：建立代表材料组织细观特征的rve模型，各组成相材料采用线弹性本构模型，在六种正交加载条件下计算其应变分布εm(x)；
[0104]
步骤s1.2：根据六种正交加载条件下的应变分布εm(x)，计算积分点处的应变集中张量a(x)；其中应变集中张量a(x)是联系宏观均质化应变εm和微观应变分布εm(x)的四阶张量，通过以下表达式计算：
[0105]
εm(x)＝a(x):εm。
[0106]
步骤s1.3：基于积分点的应变集中张量a(x)，通过k-means聚类算法对rve中的积分点划分集群，即完成以下关于应变集中张量的最小二乘优化过程：
[0107][0108]
其中，s代表最后的集群划分结果，sj代表第j个集群，an代表第n个积分点的应变集中张量，代表第j个集群内所有积分点的平均应变集中张量，n为总共的集群数量，||a||代表任意二阶张量a的2范数，即a的各个元素平方和的1/2次方。
[0109]
步骤s1.4：结合格林函数和傅里叶变换方法，计算各个集群之间的相互作用张量分量和其表达式分别为：
[0110]
[0111][0112]
其中，ω为rve的体积，ci为第i个集群的体积分数，χi(x)为第i个集群的特征函数，φ1和φ2为格林函数中与参考材料本构无关的部分，在傅里叶空间中具有以下形式：
[0113][0114][0115]
其中，ξ为傅里叶空间中的坐标向量，ξi是傅里叶空间中第i个方向的坐标，|ξ|为傅里叶空间坐标距离原点的距离，δ
ij
为kronecker delta符号，满足
[0116][0117]
参照图2所示，线上阶段包括：
[0118]
步骤s2.1：输入初始参考材料lame常数λ0和μ0，以及线下阶段生成的相互作用张量分量；
[0119]
步骤s2.2：根据参考材料lame常数λ0和μ0以及线下阶段生成的相互作用张量分量和通过以下表达式计算各个集群之间的相互作用张量d
ij
[0120][0121]
步骤s2.3：输入由离散的应力应变数据点组成的材料初始数据集，以及材料初始状态变量，同时定义材料数据集与状态变量的关系，通过数据驱动的方法求解lippmann-schwinger方程。如图3所示，该过程具体包括以下步骤：
[0122]
步骤s2.3.1：根据当前增量步的宏观积分点处的约束条件，输入微观尺度rve的均质化应变增量
[0123]
步骤s2.3.2：根据初始材料数据集e0、初始材料状态变量v0，以及t时刻第i个集群的材料状态变量v
ti
，更新t时刻第i个集群的材料数据集
[0124]
步骤s2.3.3：从材料数据集中为第i个集群分配初始的材料数据点
[0125]
步骤s2.3.4：读取t时刻第i个集群的应力应变状态以集群应力应变状态与数据点应力应变状态的能量偏差最小作为优化目标，将目标函数定义为：
[0126][0127]
其中，δεi和δσi为当前增量步的集群应变增量和集群应力增量；n为总共的集群
数量，ωi为第i个集群的体积分数，ci为第i个集群的自定义刚度矩阵。
[0128]
对应的约束条件为lippmann-schwinger方程
[0129][0130]
及均质化应变条件
[0131][0132]
其中，δε0当前增量步的参考材料应变增量；δεi和δσi为当前增量步的集群应变增量和集群应力增量；c0表示参考材料的刚度矩阵。
[0133]
步骤s2.3.5：利用拉格朗日乘子法求解上述约束条件下的目标函数最小化问题，包含约束条件的拉格朗日函数表示为：
[0134][0135]
其中，ηi和η
n+1
为拉格朗日乘子。在拉格朗日函数l中，优化变量包括δεi、δσi、δε0、ηi和η
n+1
五部分。令拉格朗日函数l对各优化变量的偏导数为零，可得到关于各优化变量的线性方程组：
[0136][0137][0138][0139][0140][0141]
其中，δ
ij
表示kronecker delta符号；ωj表示第j个集群的体积分数；cj表示第j个集群的自定义刚度矩阵；u4表示四阶单位张量；ηj表示第j个集群的拉格朗日乘子。
[0142]
步骤s2.3.6：求解以上线性方程组，获得满足约束条件的集群应变增量δεi和集群应力增量δσi；则t+1时刻集群的平衡约束点可表示为：
[0143]
[0144][0145]
实现从集群材料数据点到集群平衡约束点的映射。
[0146]
步骤s2.3.7：根据集群平衡约束点在t时刻集群材料数据集中找到能量偏差最近的数据点，集群的材料数据点可表示为：
[0147][0148]
实现从集群平衡约束点到集群材料数据点的映射。
[0149]
步骤s2.3.8：重复步骤s2.3.6～步骤s2.3.7，直到集群材料数据点达到收敛，得到t+1时刻集群平衡约束点和t+1时刻集群材料数据点
[0150]
步骤s2.3.9：根据t+1时刻集群材料数据点更新t+1时刻的集群材料状态变量v
t+1
；
[0151]
步骤s2.3.10：根据t+1时刻集群平衡约束点更新t+1时刻rve均质化应力
[0152]
步骤s2.4：根据t+1时刻集群集群平衡约束点计算rve均质化lame常数λ
rve
和μ
rve
。
[0153]
步骤s2.5：根据rve均质化lame常数(λ
rve
,μ
rve
)，更新参考材料lame常数(λ0,μ0)；重复步骤s2.2～s2.4，直到(λ
rve
,μ
rve
)与(λ0,μ0)保持一致，得到t+1时刻最终集群应力应变分布和rve均质化应力
[0154]
步骤s2.6：将rve均质化应力返回宏观积分点，获得材料的宏微观力学响应，实现加载历史相关的数据驱动多尺度并发模拟。
[0155]
下面以单向纤维增强复合材料为例，通过轴向拉伸-压缩-拉伸循环加载过程的多尺度并发模拟，对该方法进行验证。
[0156]
单向纤维增强复合材料由基体相和纤维相组成。基于matlab平台编写相关程序，通过基体材料和纤维材料的唯象本构模型，分别生成大量应力应变数据点。其中，基体材料的数据集通过弹塑性本构模型生成，弹性阶段采用各向同性线弹性模型，塑性阶段采用mises屈服准则及等向线性硬化模型。x方向的应变数据范围设置为(-0.2,0.2)，其它方向的应变数据范围设置为(-0.01,0.01)。纤维材料的数据集通过正交各向异性线弹性本构模型生成，x方向的应变数据范围设置为(-0.1,0.1)，其它方向的应变数据范围设置为(-0.01,0.01)。对于基体相和纤维相，在各自的应变数据范围内各生成10万个数据点。
[0157]
基于本方法编写子程序，通过abaqus有限元平台实现复合材料循环加载过程的数据驱动多尺度并发模拟。多尺度并发模拟过程中的宏观有限元模型和微观rve模型如图4和图5所示。为测试方便，宏观模型只包含一个单元，单元尺寸1mm，加载条件为沿x方向的拉伸-压缩-拉伸循环加载，加载参考点的位移-时间曲线如图6所示，位移幅值为0.05mm。微观
模型为代表复合材料纤维分布特征的rve模型，纤维直径为7μm，纤维体积分数为50％，rve的尺寸为50μm。为了提高计算效率，rve总共被划分为两个集群，其中基体材料和纤维材料各自被划分为一个集群。
[0158]
循环加载最终时刻宏观有限元模型及微观rve模型中沿x方向的应力分布如图7和图8所示。可以看出，在微观rve模型中，由于纤维具有更高的弹性模量，拉伸应力主要由纤维承担，基体的应力水平很小，这与轴向拉伸条件下复合材料的一般应力分布特征相符合。在宏观有限元模型中，均质化拉伸应力近似为纤维拉伸应力的一半，这与轴向拉伸条件下均质化理论的预测结果保持一致，此结果验证了本方法的准确性。
[0159]
循环加载过程中基体和纤维的应力-应变曲线如图9和图10所示。可以看出，对于纤维材料，其变形过程始终为弹性变形，其应力应变关系保持为线性关系，与加载历史无关。而对于基体材料，在变形过程中包含了塑性变形和硬化过程，其应力应变关系明显依赖于其加载历史。在第一次拉伸过程中，基体材料首先沿线性路径发生弹性变形，之后进入正向屈服和塑性硬化阶段。在第一次压缩过程中，基体材料首先沿着弹性路径卸载，之后进入反向加载过程，反向加载过程同样经历了弹性变形阶段和塑性硬化阶段，受第一次拉伸过程的影响，反向加载的初始屈服应力也有所提升。在第二次拉伸过程中，由于经历了两次硬化过程，其初始屈服应力要显著高于第一次拉伸过程。此结果证明了本发明方法可以准确预测不同尺度下材料加载历史相关的力学行为。
[0160]
本发明实施例提供了一种加载历史相关的数据驱动多尺度并发模拟方法及系统，能够在无本构模型的条件下实现材料的多尺度并发模拟，从根本上避免多尺度分析中材料本构模型所导致的系统偏差，同时能够准确预测不同尺度下材料加载历史相关的力学行为。
[0161]
本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统及其各个装置、模块、单元以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统及其各个装置、模块、单元以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同功能。所以，本发明提供的系统及其各项装置、模块、单元可以被认为是一种硬件部件，而对其内包括的用于实现各种功能的装置、模块、单元也可以视为硬件部件内的结构；也可以将用于实现各种功能的装置、模块、单元视为既可以是实现方法的软件模块又可以是硬件部件内的结构。
[0162]
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本技术的实施例和实施例中的特征可以任意相互组合。