复合材料多尺度并发模拟方法和系统

本发明涉及复合材料设计技术领域，具体地，涉及一种复合材料多尺度并发模拟方法和系统，尤其涉及一种基于自适应聚类分析及数据驱动算法的复合材料多尺度并发模拟方法

背景技术：

随着计算科学的高速发展及计算机性能的不断提升，有限元法等数值模拟技术越来越多的应用到各类工程技术领域，以指导工业生产，缩短设计周期，降低研发成本。目前流行的cae分析软件种类繁多，建模技术与求解算法日趋成熟，对于实际工程中的结构分析、高速冲击、流固耦合、多物理场等问题提供了丰富的集成模块和解决方案。而另一方面，对于大部分cae分析软件而言，其内置的材料模型不够丰富和完善，仍然存在多种形式的缺陷，特别是对于一些力学行为复杂的非线性材料，利用现有的材料模型无法准确描述其复杂的非线性力学行为。因此，准确而高效的材料本构模型，是在科学及工程问题中制约仿真分析精度及求解效率最为关键的一个因素。

复合材料作为一种典型的多相异质材料，其细观结构对宏观性能具有显著的影响。例如对于纤维增强复合材料而言，纤维的朝向、分布及体积含量，基体的塑性、损伤等非线性力学行为，以及界面的脱粘等因素都会对复合材料的宏观力学性能造成很大影响。传统的唯象本构模型从宏观尺度对材料的力学行为进行整体表征，无法描述细观尺度各相成分的应力应变演化过程；另一方面，复合材料复杂的细观结构决定其在宏观尺度往往具有明显的各向异性及强烈的非线性等特点，构造一个能够精确描述其力学行为的唯象本构模型十分困难。

代表性体积单元(representativevolumeelement,rve)方法是解决以上困难的一种有效手段，通过建立一个代表复合材料局部细观特征的周期性有限元模型来预测材料在细观尺度上力学响应演化过程。但是在多尺度并发计算时，需要在宏观有限元模型的每一个积分点处调用rve参与计算，现有的rve计算效率无法满足实际需求。自适应聚类分析方法是由美国西北大学的wingkamliu课题组提出的一种针对rve的降阶方法，将rve划分成若干集群，同一个集群具有均匀的应力和应变分布，以此降低系统的自由度数量。该方法能够极大的提升计算效率，缩短计算时间，满足多尺度并发计算需求。然而在细观尺度上，该方法仍然需要定义各个组成相的材料模型，采用唯象模型所带来的系统偏差仍然存在，在多尺度分析过程中偏差甚至可能会进一步放大，影响最终求解精度，该方法目前还难以做到完全实现数据驱动。

数据驱动算法是由加州理工学院的ortiz等人提出的一种新型计算模式，在近年来不断得以发展完善，相关理论体系已扩展至粘弹性、弹塑性及动力学等问题。该方法直接从实验所获得的材料数据集出发，结合相容性条件及平衡方程等约束条件，通过对所有材料点自由能最小化，从而获得体系的最优解。该方法完全绕过了传统计算所需的材料本构建模步骤，能够从根本上解决唯象本构模型所带来的人为误差和复杂材料性能难以建模预测的困难。然而另一方面，该方法作为一种新型计算范式，其理论体系尚不及有限元法成熟，有较多的理论问题亟需解决完善；特别是该方法目前还无法考虑复合材料多尺度并发的模拟预测，只能进行宏观结构的力学性能仿真模拟；同时由于该方法的理论体系独立于现有的有限元理论体系，因此该方法难以在现有的商用有限元软件平台上实现，基于该方法的求解器也需要独立开发、调试及维护，这给该方法的推广和应用造成了很大困难。

经过检索，专利文献cn106066913b公开了一种复杂复合材料结构等效材料性能多尺度计算方法，采用尺度分离的方法，将宏观、细观、微观三尺度结构分离，根据不同尺度模型的几何特征，分别建立各个尺度分析模型；将三尺度问题转化为两个多尺度问题：宏观-细观多尺度问题、细观-微观多尺度问题，依次对着两个多尺度问题进行分析，将微观多尺度问题得到的等效模量最终返回给宏观多尺度问题。该现有技术的不足之处在于仍然需要定义不同尺度模型，并且计算量大的技术问题。

因此，亟需研发设计一种对复合材料进行多尺度分析时，不需要引入任何经验性的材料模型的方法和系统。

技术实现要素：

针对现有技术中的缺陷，本发明的目的是提供一种复合材料多尺度并发模拟方法和系统，利用数据驱动算法构建自适应聚类分析中的各个组成相的材料本构，以数据代替模型，以便在现有的商用有限元软件平台上实现复合材料的多尺度并发计算，同时避免传统本构模型在多尺度分析中所造成的系统偏差，解决本构建模困难的难题，提高复合材料多尺度分析预测精度及计算效率。

根据本发明提供的一种复合材料多尺度并发模拟方法，包括如下步骤：

步骤s1：通过数据驱动算法构建复合材料各个组成相材料的本构；

步骤s2：针对本构模型进行自适应聚类分析，通过自适应聚类分析实现复合材料的多尺度并发计算，自适应聚类分析包括线上阶段和线下阶段，通过线下阶段基于应变集中张量对rve划分集群并求解集群之间的相互作用张量分量，通过线上阶段自适应迭代调整参考材料刚度，结合各组成相材料本构模型，获得集群应力应变分布。

优选地，步骤s1包括如下步骤：

步骤s1.1：获取复合材料各个组成相在不同加载条件下的应力应变数据集；

步骤s1.2：对于每种组成相材料，在相空间内定义代表材料自由能的罚函数形式，并确定相应的约束条件；

步骤s1.3：通过拉格朗日乘子法求解在给定约束条件下罚函数的最小值，从而确定各个组成相的材料本构。

优选地，步骤s2中线下阶段需要建立复合材料的高保真细观rve模型，计算其在不同加载条件下的弹性力学响应。

优选地，线下阶段需要计算积分点的应变集中张量，通过k-means聚类算法对积分点进行集群划分。

优选地，线下阶段通过格林函数和傅里叶变换方法计算集群之间的相互作用张量分量。

优选地，线上阶段通过参考材料刚度及相互作用张量分量计算集群之间的相互作用张量。

优选地，线上阶段输入数据驱动算法中获得的各个组成相的材料本构数据，结合各个组成相的材料本构，通过newton-raphson迭代法求解离散条件的lippmann-schwinger方程，计算所有集群的应力应变分布。

优选地，根据rve的宏观刚度调整参考材料刚度，重复线上阶段的步骤，直至宏观刚度和细观刚度一致时退出循环，得到最终的集群应力应变分布。

根据本发明提供的一种复合材料多尺度并发模拟系统，采用上述的复合材料多尺度并发模拟方法对复合材料进行多尺度并发模拟。

优选地，包括：

模块m1：通过数据驱动算法构建复合材料各个组成相材料的本构；

模块m2：针对本构模型进行自适应聚类分析，通过自适应聚类分析实现复合材料的多尺度并发计算，自适应聚类分析包括线上阶段和线下阶段，通过线下阶段基于应变集中张量对rve划分集群并求解集群之间的相互作用张量分量，通过线上阶段自适应迭代调整参考材料刚度，结合各组成相材料本构模型，获得集群应力应变分布。

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

1、本发明在对复合材料进行多尺度分析时不需要引入任何经验性的材料模型，省略了繁琐的参数标定过程，应用更为方便。

2、本发明能够利用现有的商用有限元软件平台实现多尺度并发计算，无需自主开发对应的求解器，能够有效降低计算成本

3、本发明能够直接基于复合材料各个组成相的实验数据预测复合材料的宏观力学响应，避免了传统的唯象本构模型在细观尺度及宏观尺度所导致的系统偏差，具有更高的计算精度和求解效率。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1为本发明中数据驱动算法的流程图；

图2为本发明中自适应聚类分析线下阶段流程图；

图3是自适应聚类分析线上阶段流程图；

图4是纤维材料应力应变数据集；

图5是基体材料应力应变数据集；

图6是v形缺口试样有限元模型及边界条件；

图7是v形缺口试样宏观及微观尺度剪应力分布云图。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。

如图1-3所示，本发明提供了一种复合材料多尺度并发模拟方法，包括如下步骤：

步骤s1：通过数据驱动算法构建复合材料各个组成相材料的本构。在现有的数据驱动算法中，约束条件由体系平衡方程和相容性条件来确定，通过约束条件下的极值问题直接获得体系的数值解，这使得数据驱动算法框架完全独立于现有的有限元框架，难以与有限元方法及自适应聚类分析方法相结合。而在本发明中，约束条件则用来确定组成相材料的本构关系，从而使数据驱动算法能够应用于有限元方法及自适应聚类分析方法中。

其中，步骤s1具体包括如下步骤：

步骤s1.1：获取复合材料各个组成相在不同加载条件下的应力应变数据集。以其中一种组成相材料为例，设其应力应变数据集为e，所包含的数据点为{(εi,σi),i＝1,…,n}，其中n为该材料数据集中所包含数据点的个数。不同的加载条件包括拉伸、压缩、剪切等加载条件。

步骤s1.2：对于每种组成相材料，在相空间内定义代表材料自由能的罚函数形式，并确定相应的约束条件。对于该材料相空间中的任意一个参考点z＝(ε,σ)，定义代表材料自由能的罚函数的具体形式f(z)。一般来讲，随着参考点与材料数据集两者偏差的增大，罚函数的值应增大。罚函数的具体形式并不唯一，例如可将罚函数的形式取为参考点与材料数据集中所有数据点距离平方的加权平均，如(1)式所示

其中pi为相空间中数据点(εi,σi)所对应权函数，满足归一化条件

步骤s1.3：通过拉格朗日乘子法求解在给定约束条件下罚函数的最小值，从而确定各个组成相的材料本构。

首先，根据需要确定相空间的约束条件c。对于不需要内变量的问题，如线弹性问题，需要求解材料在任意给定应变状态ε0条件下的应力，则相应的约束条件为：

c＝{(σ,ε),ε＝ε0}(2)

对于需要内变量的问题，如弹塑性问题，需要求解材料在任意给定应变状态ε0和内变量状态q0条件下的应力，则相应的约束条件为：

c＝{(σ,ε,q),ε＝ε0,q＝q0}(3)

其次，求解罚函数在约束条件下的极小值问题以(1)式所示的罚函数形式和(2)式所示的约束条件为例，通过拉格朗日乘子法，该约束极值问题可转换为以下形式的驻值问题

通过求解(4)式可获得在任意应变ε0约束条件下材料在相空间的位置(σ0,ε0)，即为材料的本构关系。

最后，通过数值方法计算材料的切线刚度矩阵s。在当前的应力应变状态(σ0,ε0)条件下，分别给应变的六个分量一个小增量dε，利用以上方法计算在六种应变状态ε0+dε下的应力状态，得到对应的六个应力增量dσ，将其进行归一化组装即可得到材料的切线刚度矩阵。

步骤s2：针对本构模型进行自适应聚类分析，通过自适应聚类分析实现复合材料的多尺度并发计算。自适应聚类分析主要包括线下阶段和线上阶段两部分，其主要流程分别如图2和图3所示。线下阶段主要基于应变集中张量对rve划分集群，以此降低系统自由度数量，并求解集群之间的相互作用张量分量；线上阶段则通过自适应迭代调整参考材料刚度，结合各组成相材料本构模型，通过求解lippmann-schwinger获得集群应力应变分布。

对于线下阶段，具体如下：

(1)建立复合材料的高保真细观rve模型，高保真细观rve模型要能够反映复合材料的局部细观结构特征。通过计算rve在六种正交载荷条件下的弹性力学响应，获得弹性应变分布ε^micro(x)，之后根据表达式ε^micro(x)＝a(x):ε^macro(x)计算应变集中张量分布a(x)；

(2)基于应变集中张量通过k-means聚类算法对rve进行集群划分，使应变集中张量接近的积分点划分至同一个集群内，在同一个集群认为具有均匀的应力和应变，以此降低系统自由度数量；

(3)利用格林函数和傅里叶变换等方法，计算集群之间的相互作用张量分量和其中相互作用张量分量和为相互作用张量d^ij中不随参考材料本构变化而变化的部分。

对于线上阶段，具体如下：

(1)给定参考材料lame常数λ⁰和μ⁰，通过相互作用张量分量和计算相互作用张量d^ij，计算表达式如(5)式所示

(2)根据数据驱动算法中的本构模型，计算在给定第j个集群的应变增量δε^j的条件下，其应力增量δσ^j和切线刚度矩阵c^j；

(3)将应力增量δσ^j和切线刚度矩阵c^j代入离散条件下的lippmann-schwinger方程(6)中

通过newton-raphson迭代法求解集群的应变增量δε^j，更新集群应力应变分布状态；

(4)根据集群应力应变增量δσ^j和δε^j，计算rve的宏观刚度矩阵所对应的lame常数λ和μ，并与参考材料lame常数λ⁰和μ⁰相对比。若两者相对偏差超出给定范围，则将(5)式中参考材料lame常数λ⁰和μ⁰更新为λ和μ，并重复步骤(1)-(3)，直至两者相对偏差在给定范围之内时退出循环，此时集群应力应变分布即为所求rve的应力应变分布状态。

下面以纤维增强复合材料为例，通过v形缺口试样剪切测试对该方法进行验证。

纤维增强复合材料主要由纤维增强相与树脂基体相组成，纤维和基体的应力应变数据集分别如图4和图5所示；v形缺口试样有限元模型及加载条件如图6所示，左夹持端固定，右夹持端施加向上的位移载荷1mm。

基于本发明编写了vumat子程序，通过abaqus有限元平台实现了v形缺口试样在剪切条件下的多尺度并发计算。v形缺口试样最终时刻的宏观及微观剪应力分布云图如图7所示，可以看出在宏观尺度上v形缺口试样的高应力区主要集中在中心区域，在微观尺度上剪应力主要由纤维材料承担，基体材料承担的剪应力相对较小。该算例验证了本专利方法的有效性。

本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统及其各个装置、模块、单元以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统及其各个装置、模块、单元以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同功能。所以，本发明提供的系统及其各项装置、模块、单元可以被认为是一种硬件部件，而对其内包括的用于实现各种功能的装置、模块、单元也可以视为硬件部件内的结构；也可以将用于实现各种功能的装置、模块、单元视为既可以是实现方法的软件模块又可以是硬件部件内的结构。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。