CN109241562B - 基于多尺度有限元方法的微结构材料弹性性能测定方法 - Google Patents

Abstract

一种基于多尺度有限元方法的微结构材料弹性性能测定方法，根据超材料单胞结构形式与尺寸构建检测模型并对检测模型进行粗网格划分，并在划分后的网格上的突变位置附加虚拟节点；然后在得到的单胞上进行子网格划分，并建立宏观粗网格单元内部的子网格节点的位移与所有节点间的多尺度基函数关系，根据多尺度基函数推导得到粗网格单元刚度矩阵并组装得到整体刚度矩阵，从而获得宏观位移信息，即微结构材料弹性性能。本发明能够实现显著降低多尺度有限元方法应用于微结构材料误差的效果，解决单胞与单胞间连接不匹配的计算问题。

Description

基于多尺度有限元方法的微结构材料弹性性能测定方法

技术领域

本发明涉及的是一种增材制造领域的技术，具体是一种基于多尺度有限元方法的微结构材料弹性性能测定方法。

背景技术

增材制造(AM)是采用材料累加方法制造实体零件的技术，在金属与非金属材料制造领域中均得到众多成果与应用。由于增材制造可以制造复杂造型与细微结构均这一重要特点，面向增材制造的优化设计分为面向宏观复杂结构外形的宏观外形优化与面向细微结构的微结构优化两个领域。

在微结构材料优化领域中，传统多尺度有限元方法并不是针对增材制造的微结构材料开发的，且单胞间连接问题没有考虑，所以计算仍具有较大误差。针对增材制造微结构弹性性能预测问题，如何利用嵌入式协同的多尺度有限元方法，构建适合于增材制造微结构仿真的精确多尺度有限元模型，同时考虑微结构各单胞间的连接问题，是保证预测精度与效率的关键。

发明内容

本发明针对现有技术存在的上述不足，提出一种基于多尺度有限元方法的微结构材料弹性性能测定方法，在构建宏观有限元单元刚度矩阵的过程中，不同于常规多尺度有限元通过固定节点传递多尺度信息的特点，而是通过将虚拟节点与边界材料性质突变点结合，从而将微观单胞结构与宏观粗网格节点关联，达到显著降低多尺度有限元方法应用于微结构材料误差的效果。通过考虑单胞与相邻单胞的结构特征附加粗网节点，计算过程中考虑单元连接情况，从而解决了单胞与单胞间连接不匹配的计算问题。

本发明是通过以下技术方案实现的：

本发明涉及一种基于多尺度有限元方法的微结构材料弹性性能测定方法，根据超材料单胞结构形式与尺寸构建检测模型并对检测模型进行粗网格划分，并在划分后的网格上的突变位置附加虚拟节点；然后在得到的单胞上进行子网格划分，并建立宏观粗网格单元内部的子网格节点的位移与所有节点，即粗网格原有节点与虚拟节点间的多尺度基函数关系，根据多尺度基函数推导得到粗网格单元刚度矩阵并组装得到整体刚度矩阵，从而获得宏观位移信息，即微结构材料弹性性能。

所述的粗网格划分是指：根据选用单胞形状、尺寸在检测模型的设计域上确定粗网格划分形式，优选划分为四节点正方形单元，单元尺寸与所用单胞尺寸一致，划分的网格所具有的节点为固有节点。

所述的虚拟节点是指：在每个粗网格边界上对应单胞边界处材料性质突变位置附加的节点，其位置随单胞结构边界形式变化而变化，从而实现单胞结构的边界形式与虚拟节点形式一一对应。

所述的虚拟节点进一步经过重合判断，由于相邻单胞共用一条粗网格边界，边界两边的单胞均在边界上虚拟节点，当虚拟节点位置重合时，重合的两个节点合并为一个节点，否则保留全部节点。

所述的子网格划分是指：在虚拟节点和固有节点的基础上得到单胞，对每个单胞进行子网格划分，其形式与尺寸依据单胞结构特征确定，然后建立宏观粗网格单元内部子网格节点位移与粗网格个节点间的多尺度基函数关系。

所述的子网格采用但不限于三角形网格、四边形网格或四面体网格、六面体网格等。

所述的多尺度基函数关系是指：u_{e_all}＝Nu_E，其中：u_{e_all}为粗网格内部子网格上所有节点的位移向量；u_E为宏观粗网格单元节点的位移向量；N为各个节点基函数组成的矩阵：

其中：k为子网格节点数目。该基函数对于任何子网格节点，当采用普通多尺度有限元思想的表达式如下：

当采用扩展的多尺度有限元思想，基函数表达形式如下：

其中：n表示步骤二中确定的该单胞对应粗网格节点数目总和。

所述的基函数选用扩展的多尺度有限元思想构造基函数，其进一步满足粗网格的刚体位移与单元间的位移协调，即：

从而进一步考虑到不同方向的力与位移的耦合关系，在一定程度上提高计算精度。

所述的粗网格单元刚度矩阵

其中：e为粗网格中任意一个细网格；m为细网格数目；K_E为粗网格单元刚度矩阵；L_e为该粗网格单元与对应子网格单元之间的映射关系矩阵，由多尺度基函数构成：

其中下标_ej指的是第_e个子网格中的第_j个节点；D_e和B_e分别为子网格单元的应变矩阵与弹性矩阵，即常规有限元方法中使用的应变矩阵与弹性矩阵。

所述的整体刚度矩阵为：

其中P指划分的粗网格总数，K_Ei为第i个粗网格对应的刚度矩阵，K是由每个粗网格单元刚度矩阵K_Ei在同一总体坐标系下按照每个节点编号对应位置叠加而成的。

所述的宏观位移信息是指：通过整体刚度矩阵和位移边界条件与力的边界条件求解平衡方程，具体为：Kδ＝F，其中δ为粗网格节点位移向量，F为载荷向量。

所述的宏观计算结果可以用来进行微观层面后处理，从而获得微结构内部应力应变等信息，其采用但不限于：粗网格节点位移向量δ_E由平衡方程求解获得，子网格单元节点位移δ_e由粗网格单元与对应子网格单元之间的映射关系矩阵得到：δ_e＝L_eδ_E。从而可以获得该子网格单元内部应力：ζ_e＝D_eB_eδ_e；该子网格单元内部应变：ε_e＝B_eδ_e。

本发明涉及一种实现上述方法的系统，包括：粗网格划分模块、节点附加模块、子网格划分模块、基函数模块、刚度矩阵模块以及整体刚度预测模块，其中：粗网格划分模块与节点附加模块相连并传输单胞边界信息，粗网格划分模块与子网格划分模块相连并传输单胞划分信息，子网格划分模块与基函数模块相连并传输子网格节点信息，基函数模块与刚度矩阵模块相连并传输粗网格节点与细网格节点位移关联信息，刚度矩阵模块与整体刚度预测模块相连并传输粗网格刚度信息。

技术效果

与现有技术相比，本发明基于多尺度有限元方法，宏观粗网格节点依赖于单胞结构形式或材料性质突变位置，自动调整粗网格节点位置，宏观节点与微观节点关联，在布置少量节点的情况下极大提升计算精度，同时保证计算效率。传统均匀化方法与多尺度有限元方法假设单胞间连接连续，并不适合于此类微结构问题的计算，而本发明在附加粗网格节点时，不仅考虑本单胞结构，还考虑了相邻单胞结构形式，进而可以对连接不连续的问题进行准确仿真，扩展了多尺度有限元在微结构仿真计算中的应用。传统基于均匀化方法的优化问题往往只能利用密度作为优化变量，建立材料密度与材料矩阵之间的联系，但实际微结构形式多样，密度无法准确描述材料性质。采用本发明可以利用微结构多种相关几何特征作为优化变量，将微结构力学性质与粗网格单元刚度矩阵关联，建立起微结构特征与材料性质的单一联系，对于使用多种形式单胞进行准确的优化提供技术支持。

附图说明

图1为本发明方法流程示意图；

图2为本发明中微结构超材料单胞示意图；

图3为本发明中粗网格划分示意图；

图4为本发明中设计域内单胞分布示意图；

图5为本发明中设计域边界条件示意图；

图6为本发明中有限元精细网格建模示意图；

图7为本发明中结果对比图。

具体实施方式

如图1所示，本实施例以两种平面矩形超材料间隔分布构件的弯曲变形为例进行说明，具体包括以下步骤：

步骤一、根据选用单胞形状、尺寸在设计域划分粗网格，具体步骤包括：

1.1)确定本次实施例中选用的超材料单胞结构形式与尺寸，本实施例针对二维问题，选择两种单胞类型如图2(a)所示，两种单胞所对应正方形边长均为l，单胞一整体表现为正泊松比，单胞二为负泊松比。

1.2)选择设计域为5×2的矩形区域，依据选择的单胞，将设计域进行粗网格划分，划分结果如图3(a)所示，粗网格与单胞尺寸一致。粗网格划分为4节点矩形单元，这些节点即为固定节点，仅仅与单胞外形尺寸相关，与内部结构无关。

如图4所示。在优化过程中，设计域内单胞分布形式不断变动，在本实施例中选定某种尺寸的单胞，在一个设计周期内将不再变动，所以粗网格往往只需要划分一次，所以固定节点位置不变。当选区单胞尺寸发生变化，粗网格需要重新划分。

步骤二、根据子各个单胞结构形式确定粗网格虚拟节点位置，具体为：

2.1)基于步骤一中的粗网格划分，单胞外界被一个粗网格包围。沿此边界进行搜索，当遇到材料性质发生突变，记录此位置，并在响应粗网格边界虚拟节点。

由于本实施例当中组成微结构的材料设定为单一且均匀连续的，即为搜索边界上的孔洞与材料交界位置。对于两种单胞形式，搜索结果如图2(b)所示，编号为1、2、3、4的节点即为固定节点，其余节点为搜索得到的虚拟节点。

2.2)将虚拟节点放置于粗网格上。由于相邻单胞共用粗网格边界，当两边虚拟节点位置重合，将造成计算资源的浪费。由记录的虚拟节点位置，判别是否出现虚拟节点重合，若重合则将两节点合并。进行节点附加的粗网格如图3(b)所示。

步骤三、子网格划分，并计算多尺度基函数，具体为：

3.1)依据单胞的结构特点，选取合适的单元进行子网格划分，本实施例中选取四节点矩形单元进行划分，由步骤2中获取的虚拟节点位置信息，在单胞相应位置布置子网格节点，保证粗细网格间信息传递的一致性。边界处由于考虑相邻单胞虚拟节点影响，需要将于边界相邻的无材料区域同样划分网格。网格划分效果如图2(c)所示。选择单胞材料进行材料矩阵计算，选择钢材作为研究对象，将深色区域网格赋予钢材力学性能参数，浅色网格赋予极小的弹性模量参数来模拟空材料区域。子网格单元划分信息与材料力学性质基本参数如表1和表2所示。

表1子网格单元划分信息

表2材料力学性能参数

3.2)利用获得的子网格进行多尺度基函数计算。依据设计要求、试件情况选择基函数构造形式与边界条件形式，采用扩展的多尺度有限元思想构造多尺度基函数，选择线性边界条件进行计算。对于子网格具有m个节点的单胞i，由步骤一中得到4个固定节点与步骤二中搜索得到的n_i个虚拟节点，通过附加MsFEM所需的边界条件求解平衡方程，对于此二维问题相当于求解(4+n_i)×2次平衡方程，得到子网格中任意节点与粗网格上节点之间的关系，即基函数。该单胞基函数数据量为2m(4+n_i)个。

步骤四、计算粗网格的等效刚度矩阵，并根据边界条件进行宏观计算，具体为：

4.1)由步骤三得到的全部节点基函数，挑选出针对各个节点的基函数，并组成节点基函数矩阵L_e，应变矩阵B_e与网格划分相关，材料矩阵D_e与该子网格材料相关；利用能量守恒原理即可求得宏观粗网格单元的刚度矩阵，此刚度矩阵包含了单胞内部结构信息(即基函数)，同时由于此刚度阵附加了相邻单胞的节点，可以计算单胞间不连续的情况。

所述的刚度矩阵

4.2)对每一个粗网格的刚度矩阵进行计算，并按照节点编号组装成整体刚度矩阵。本问题边界条件如图5所示，设计域左侧限制x、y方向位移，右侧受到均匀分布向上的力，大小为96。由整体刚度矩阵与引入的边界条件，求解平衡方程，得到设计零件的位移信息。

结果验证：将普通有限元精细建模方法求得该问题的解与通过以上四步所得到微结构材料宏观力学响应比较，普通有限元方法网格划分如图6所示。同时将均匀化方法、传统的多尺度有限元、传统的多尺度有限元法、扩展的多尺度有限元方法对该问题的计算结果与本方法进行对比，对比如图7所示。对比图展示的是零件y挠度图，可以看出，运用本方法与有限元精确建模方法基本一致，而传统多尺度方法对于这类问题计算误差较大。本方法与两种传统多尺度方法对于零件最右侧挠度误差对比如表3所示。

表3三种方法与精细建模误差对比

上述具体实施可由本领域技术人员在不背离本发明原理和宗旨的前提下以不同的方式对其进行局部调整，本发明的保护范围以权利要求书为准且不由上述具体实施所限，在其范围内的各个实现方案均受本发明之约束。

Claims (10)

1.一种基于多尺度有限元方法的微结构材料弹性性能测定方法，其特征在于，根据超材料单胞结构形式与尺寸构建检测模型并对检测模型进行粗网格划分，并在划分后的网格上的突变位置附加虚拟节点；然后在得到的单胞上进行子网格划分，并建立宏观粗网格单元内部的子网格节点的位移与所有节点间的多尺度基函数关系，根据多尺度基函数推导得到粗网格单元刚度矩阵并组装得到整体刚度矩阵，从而获得宏观位移信息，即微结构材料弹性性能；

所述的虚拟节点是指：在每个粗网格边界上对应单胞边界处材料性质突变位置附加的节点，其位置随单胞结构边界形式变化而变化，从而实现单胞结构的边界形式与虚拟节点形式一一对应；

所述的多尺度基函数关系是指：u_{e_all}＝Nu_E，其中：u_{e_all}为粗网格内部子网格上所有节点的位移向量；u_E为宏观粗网格单元节点的位移向量；N为各个节点基函数组成的矩阵：

其中：k为子网格节点数目；基函数对于任何子网格节点，当采用普通多尺度有限元思想的表达式如下：

当采用扩展的多尺度有限元思想，基函数表达形式如下：

其中：n表示单胞对应粗网格节点数目总和。

2.根据权利要求1所述的方法，其特征是，所述的粗网格划分是指：根据选用单胞形状、尺寸在检测模型的设计域上确定粗网格划分形式，单元尺寸与所用单胞尺寸一致，划分的网格所具有的节点为固有节点。

3.根据权利要求1所述的方法，其特征是，所述的虚拟节点进一步经过重合判断：当虚拟节点位置重合时，重合的两个节点合并为一个节点，否则保留全部节点。

4.根据权利要求1所述的方法，其特征是，所述的子网格划分是指：在虚拟节点和固有节点的基础上得到单胞，对每个单胞进行子网格划分，其形式与尺寸依据单胞结构特征确定，然后建立宏观粗网格单元内部子网格节点位移与所有节点，即粗网格原有节点与虚拟节点间的多尺度基函数关系。

5.根据权利要求1所述的方法，其特征是，所述的采用扩展的多尺度有限元思想构造基函数，其进一步满足粗网格的刚体位移与单元间的位移协调，即：

从而进一步考虑到不同方向的力与位移的耦合关系并提高计算精度。

6.根据权利要求1所述的方法，其特征是，所述的粗网格单元刚度矩阵

其中：e为粗网格中任意一个细网格；m为细网格数目；K_E为粗网格单元刚度矩阵；L_e为该粗网格单元与对应子网格单元之间的映射关系矩阵，由多尺度基函数构成：

其中下标ej指的是第e个子网格中的第j个节点；D_e和B_e分别为子网格单元的应变矩阵与弹性矩阵，即常规有限元方法中使用的应变矩阵与弹性矩阵。

7.根据权利要求6所述的方法，其特征是，所述的整体刚度矩阵为：

其中P指划分的粗网格总数，K_Ei为第i个粗网格对应的刚度矩阵，K是由每个粗网格单元刚度矩阵K_Ei在同一总体坐标系下按照每个节点编号对应位置叠加而成的。

8.根据权利要求1所述的方法，其特征是，所述的宏观位移信息是指：通过整体刚度矩阵和位移边界条件与力的边界条件求解平衡方程，具体为：Kδ＝F，其中δ为粗网格节点位移向量，F为载荷向量。

9.根据权利要求1所述的方法，其特征是，所述的宏观计算结果进一步用来进行微观层面后处理，从而获得微结构内部应力应变等信息，具体为：粗网格节点位移向量δ_E由平衡方程求解获得，子网格单元节点位移δ_e由粗网格单元与对应子网格单元之间的映射关系矩阵得到：δ_e＝L_eδ_E，从而获得该子网格单元内部应力：σ_e＝D_eB_eδ_e；该子网格单元内部应变：ε_e＝B_eδ_e。

10.一种实现上述任一权利要求所述方法的系统，其特征在于，包括：粗网格划分模块、节点附加模块、子网格划分模块、基函数模块、刚度矩阵模块以及整体刚度预测模块，其中：粗网格划分模块与节点附加模块相连并传输单胞边界信息，粗网格划分模块与子网格划分模块相连并传输单胞划分信息，子网格划分模块与基函数模块相连并传输子网格节点信息，基函数模块与刚度矩阵模块相连并传输粗网格节点与细网格节点位移关联信息，刚度矩阵模块与整体刚度预测模块相连并传输粗网格刚度信息。

Images