CN113283147A - 三维Cohesive单元全局嵌入方法 - Google Patents

Abstract

本发明公开了一种三维Cohesive单元全局嵌入方法，包括以下步骤：S1，建立裂纹三维的有限元模型，网格剖分并导入实体，保存为模型文件输出；S2，读取并提取模型文件中的单元数据及节点数据，保存为数据文件输出；S3，读取数据文件，获得单元矩阵和节点矩阵；分别对单元矩阵和节点矩阵进行处理，得到扩展节点矩阵和重构单元矩阵；原单元矩阵和重构单元矩阵的节点编号列组合，得到原面矩阵和重构面矩阵；S4根据原面矩阵和重构面矩阵获得三维Cohesive矩阵；将扩展节点矩阵，重构单元矩阵以及Cohesive矩阵替换原节点数据和单元数据，得到新模型文件；S5，工程处理软件读取新模型文件，得到全局嵌入三维Cohesive单元的模型。

Description

三维Cohesive单元全局嵌入方法

技术领域

本发明涉及裂纹有限元离散方法建模分析领域，具体涉及一种三维Cohesive单元全局嵌入方法。

背景技术

随着监测技术和模拟技术的发展，学者们对于裂纹形态、裂纹扩展机理等研究已不再局限于二维平面内的欧式线性裂纹，同时，越来越复杂的几何模型和数学模型使得数值模拟等技术对于当前的研究分析工作显得越发重要。数值模拟技术能解决理论推导的解析解假定太多而不容易符合实际工程的问题。其在多尺度研究中的作用尤为明显，不仅能用于大尺度的复杂工程中局部破坏、局部应力不均等现象的描述，便于工程师们对工况和研究结果进行对比分析，在实验室研究中，模拟求解得到的历程云图、随时间变化的位移云图、孔隙压力云图等为学者们对本质现象的描述和成因分析提供了有力的参考和引导。

数值模拟技术发展至今，已经开始由有限元逐渐向离散元发展，然而，即使在计算机快速发展，且普及程度已经如此之高的时代，有限元作为一种最早发展出来的数值模拟方法，仍然具有其独特的优势和巨大的潜力。这是由于离散元等数值模拟技术的基本理论还不够充分，同时，其计算工程过于缓慢，即使是在当今，其对计算机性能的要求也略高。也就有了过渡型的一些方法，如边界元、有限离散元法等，这些模拟方法基于有限元的基本理论，同时采用离散的基本思想将地质模型、工程实体模型等抽象为由许多块体，许多节理裂隙分割的，却又保持为整体的模型，既保证了求解速度，又能够在一定程度上再现工程实体的局部响应。ABAQUS作为大型非线性有限元求解软件，一直受到极大的关注，其内裂纹的模拟手段之一的Cohesive单元建模方法与有限离散元法的思想基本一致，因而在航空航天、机械制造、材料科学、岩土工程等领域具有相当的潜力。然而，当前的文献等资料中尚未有公开提到的ABAQUS三维Cohesive单元的全局嵌入方法，使得ABAQUS三维Cohesive单元的嵌入繁琐，影响建模仿真效率，不便于对裂纹扩展数值模拟等进行研究。

发明内容

本发明的目的在于克服现有技术中所存在的不能快速全局嵌入二维Cohesive单元的不足，提供一种三维Cohesive单元全局嵌入方法，在短时间内对ABAQUS中裂缝图像的裂纹几何模型实现全局嵌入Cohesive单元，极大提高了ABAQUS的建模效率，便于后续的裂缝扩展数值模拟研究。

为了实现上述发明目的，本发明提供了以下技术方案：

一种三维Cohesive单元全局嵌入方法，包括以下步骤：

S1，采集裂纹信息，通过工程处理软件建立三维的有限元模型，进行网格剖分，并导入实体，处理后的文件保存为模型文件，并输出该模型文件；

S2，通过数据处理软件读取三维的模型文件，提取模型文件中的单元数据及节点数据，并将提取的数据保存为数据文件，并输出该数据文件；

S3，采用数据分析软件读取数据文件记作原始矩阵，原始矩阵通过矩阵分割获得单元矩阵和节点矩阵；分别对单元矩阵和节点矩阵进行处理，将各立体单元共用的节点及面区分开，分别得到扩展节点矩阵和重构单元矩阵；单元矩阵和重构单元矩阵的节点编号列组合，得到原面矩阵和重构面矩阵；

S4，根据原面矩阵获取各六面体单元共用的面，再根据重构面矩阵获得三维Cohesive矩阵；将扩展节点矩阵，重构单元矩阵以及三维Cohesive矩阵放入模型文件的对应位置替换原节点数据和单元数据，得到新模型文件；

S5，工程处理软件读取新模型文件，得到全局嵌入三维Cohesive单元的模型，进行模拟分析。

优选地，所述步骤S3包括以下步骤：

步骤S31，采用数据分析软件读取数据文件记作原始矩阵，区分原始矩阵的单元数据和节点数据，获得单元矩阵和节点矩阵；

步骤S32，对节点矩阵进行扩展，并对单元矩阵进行重新构建；对单元矩阵内的重复节点重新编号，生成重构单元矩阵，并且扩展节点矩阵包含重新编号的节点；根据原单元矩阵获取表示单元组成面的原面矩阵，根据重构单元矩阵获得表示重构单元组成面的重构面矩阵。

优选地，所述步骤S31采用MATLAB软件作为数据分析软件，首先判定原始矩阵中的非数值元素，排除非数值元素的干扰，通过矩阵分割生成节点矩阵和单元矩阵。

优选地，所述步骤S31采用isnan函数和sum函数进行矩阵分割，获得单元矩阵和节点矩阵。

优选地，所述步骤S32包括以下步骤：

步骤S321，对节点矩阵进行扩展，节点矩阵的扩展方式为，提取节点矩阵的第一列数据节点编号，将原节点编号扩展为多个不重复的扩展节点编号，使每个单元的节点独立存在，且扩展后的扩展节点对应的平面坐标值为原节点的平面坐标值；

步骤S322，对单元矩阵进行重新构建，单元矩阵的重构方式为，对单元矩阵包含的重复节点进行重新编号，使单元矩阵包含的节点不重复；重复节点的编号方式与节点矩阵的节点编号扩展方式相对应；

S323，单元矩阵和重构单元矩阵的节点编号列组合，得到原面矩阵和重构面矩阵；所述原面矩阵和重构面矩阵用于表示组成六面体的六个面。

优选地，所述步骤S322包括以下步骤：

步骤S3221，更改单元矩阵的单元编号列，区分单元编号与节点编号的数据区间；将更改单元编号列后的单元矩阵变形，生成单行或单列的单元矩阵；

步骤S3222，获取节点编号不重复的节点数目值，建立for循环，根据for循环，查找各节点在单行或单列的单元矩阵中出现的次数及位置，对单行或单列的单元矩阵包含的重复节点进行重新编号，重复节点的编号方式与节点矩阵的节点编号扩展方式相对应；单行或单列的单元矩阵的重复节点重新编号后，通过矩阵分割或变形，转换为与原始单元矩阵大小一致的重构单元矩阵。

优选地，所述步骤S3222中采用unique函数和length函数获取节点编号不重复的节点数目值Unique_C3D8L；采用unique函数去掉已删除单元编号列的单元矩阵中的重复元素，得到单值矩阵Unique_C3D8，采用length函数计算矩阵Unique_C3D8的长度Unique_C3D8L；所述步骤S3222中for循环内通过ismember函数及bwlabel函数，得到重复节点出现次数及位置；ismember函数用于将单行或单列的单元矩阵二值化，以判断单元矩阵中的数据是否为待判断的节点；bwlabel函数用于对二值化后的单元矩阵进行连通域判断，以得到重复节点出现次数及位置。

优选地，所述步骤S4包括以下步骤：

S41，根据原面矩阵查找各六面体单元所共有的面，并查找共用的面在原面矩阵和重构面矩阵中的位置，然后查询在原面矩阵和重构面矩阵中用于表示共用的面的节点编号，查询的原画矩阵的节点编号根据工程处理软件的编号规则进行重新排序；查询的重构面矩阵的节点编号用于表示三维Cohesive单元的组成节点，查询的重构面矩阵的节点编号根据原画矩阵中节点的重新排序的规则排序后生成用于添加三维Cohesive单元的三维Cohesive矩阵；

S42，删除扩展节点矩阵，重构单元矩阵以及三Cohesive矩阵多余的数据部分后，放入模型文件的对应位置替换原节点数据和单元数据，得到新模型文件。

优选地，所述步骤S41询共用面的位置的方式为：将原面矩阵中表示六面体单元各面的四列节点编号数字分别相加和相乘，之后再将获得的和与乘积以字符串进行组合，组合后将其重新转化为数字，定义为变量Product_Sum；利用unique函数查找Product_Sum中的单值，建立for循环，通过find函数寻找Product_Sum矩阵中不同变量的位置矩阵Sit_find，若位置矩阵长度为2，则其对应的单元共面，共用面为位置矩阵指示的面。

优选地，所述三维Cohesive单元为流固耦合单元或非流固耦合单元。

与现有技术相比，本发明的有益效果：根据ABAQUS输出的INP文件特点，利用EXCEL和MATLAB矩阵实验室，对ABAQUS的inp文件中的坐标和节点编号矩阵和六面体单元的节点和单元编号矩阵进行数据处理，建立全局三维Cohesive单元的节点和编号矩阵，再将节点矩阵和单元矩阵按照ABAQUS的inp文件格式修改为在ABAQUS中能够读取的包含节点和单元数据的inp文件。本发明主要利用MATLAB中的各种函数对ABAQUS中inp内的节点和单元矩阵进行数据处理，能够建立全局嵌入三维Cohesive单元的代码库，在短时间将ABAQUS中地质体几何模型实现全局嵌入三维Cohesive单元，极大提高了ABAQUS的建模速率，且采用外部程序对INP文件的几何数据进行处理，降低了ABAQUS版本变化对全局嵌入三维Cohesive单元的影响。

附图说明：

图1为本发明示例性实施例1的三维Cohesive单元全局嵌入方法的流程图一；

图2为本发明示例性实施例1的三维Cohesive单元全局嵌入方法的流程图二；

图3为本发明示例性实施例1的的ABAQUS裂纹的示例三维几何模型；

图4为本发明示例性实施例1的全局嵌入三维Cohesive单元后的Cohesive单元集合示意图；

图5为本发明示例性实施例1的六面体单元的示意图；

图6为本发明示例性实施例1的六面体单元共面的示意图；

图7为本发明示例性实施例1的三维Cohesive单元的获取流程示意图；

图8为本发明示例性实施例1的三维非流固耦合的Cohesive单元的添加示意图；

图9为本发明示例性实施例1的三维流固耦合的Cohesive单元的添加示意图。

具体实施方式

下面结合试验例及具体实施方式对本发明作进一步的详细描述。但不应将此理解为本发明上述主题的范围仅限于以下的实施例，凡基于本发明内容所实现的技术均属于本发明的范围。

实施例1

如图1或图2所示，本实施例提供一种三维Cohesive单元全局嵌入方法，根据ABAQUS中对三维Cohesive单元的定义规则和inp文件中节点和单元的排列格式，通过EXCEL读取ABAQUS输出的inp文件，提取inp文件中与节点和单元数据相关的数字矩阵，导入到MATLAB中进行矩阵的数据处理，重新建立inp文件，实现三维Cohesive单元的全局嵌入。具体包括以下步骤：

S1，采集裂纹信息，通过工程处理软件建立三维的有限元模型，进行网格剖分，并导入实体，处理后的文件保存为模型文件，并输出该模型文件；

S2，通过数据处理软件读取三维的模型文件，提取模型文件中的单元数据及节点数据，并将提取的数据保存为数据文件，并输出该数据文件；

S3，采用数据分析软件读取数据文件记作原始矩阵，原始矩阵通过矩阵分割获得单元矩阵和节点矩阵；分别对单元矩阵和节点矩阵进行处理，将各立体单元共用的节点及面区分开，分别得到扩展节点矩阵和重构单元矩阵；原单元矩阵和重构单元矩阵的节点编号列组合，得到原面矩阵和重构面矩阵；

S4，根据原面矩阵获取各六面体单元共用的面，再根据重构面矩阵获得三维Cohesive矩阵；将扩展节点矩阵，重构单元矩阵以及三维Cohesive矩阵放入模型文件的对应位置替换原节点数据和单元数据，得到新模型文件；

S5，工程处理软件读取新模型文件，得到全局嵌入三维Cohesive单元的模型，进行模拟分析。

其中，步骤S1的详细过程如下所示：

采用ABAQUS软件作为工程处理软件。如图2所示，通过ABAQUS建立三维的有限元模型，进行网格剖分，并导入模型实体，处理后的文件保存为后缀名为inp的模型文件，并输出该模型文件，所述模型文件记作C3D8_Coh.inp。

其中，步骤S2的详细过程如下所示：

采用EXCEL软件作为数据处理软件。通过EXCEL读取该INP文件(即文件C3D8_Coh.inp)，提取模型文件中的节点数据和单元数据，另存为后缀名为xlsx的数据文件，所述数据文件记作C3D8_Coh.xlsx。

其中，步骤S3的详细过程如下所示：

步骤S31，采用数据分析软件读取数据文件记作原始矩阵，区分原始矩阵的单元数据和节点数据，获得单元矩阵和节点矩阵；

具体地，数据分析软件可以为应用MATLAB、java、VB、C语言、python等多种语言的汇编软件。本实施例采用MATLAB软件作为数据分析软件。MATLAB读取步骤S2所述的xlsx文件(C3D8_Coh.xlsx)，通过矩阵分割获得节点矩阵和单元矩阵，所述节点矩阵记作Node，所述单元矩阵记作C3D8。所述节点数据的列数为四列，其中第一列为节点编号，第二列至第四列分别为节点的三维坐标；所述单元数据的列数为九列，其中，第一列是单元编号，后面八列是单元包含的八个节点的编号。

但是MATLAB读取的单元数据可能存在列错位的现象，因此需排除列错位带来的非数值元素的干扰。

优选地，首先判定原始矩阵中的非数值元素，排除非数值元素的干扰，通过矩阵分割生成节点矩阵和单元矩阵；

排除非数值元素干扰的步骤如下所示：删除每一行最后一个数值元素前的非数值元素，该行剩余元素位置前移，再通过判断处理后的矩阵中数值元素的列数，进行矩阵分割，得到节点矩阵和单元矩阵；或者识别原始矩阵中的数值块，通过判断数值块的列数，识别数值块的类别(属于节点矩阵或单元矩阵)，再将同一类别的数值块分割拼接，得到节点矩阵和单元矩阵；所述数值块中每行的列数相同且列数的序号一致。

优选地，采用isnan和sum函数获得节点矩阵和单元矩阵的行数和列数，通过矩阵分割成节点矩阵Node和单元矩阵C3D8。

数据文件中的空格输入MATLAB则表示为非数值NAN。MATLAB读取数据文件得到的矩阵后采用isnan函数判定非数值元素，排除非数值元素的干扰。通过isnan函数判定矩阵中的元素是否为无穷大。isnan函数返回一个与表格相同维数的数组，若元素为非数值NaN，在对应位置上返回逻辑1(真)，否则返回逻辑0(假)；然后排除非数值元素的影响，将读取数据文件后得到的矩阵通过矩阵分割生成节点矩阵和单元矩阵。排除非数值元素干扰的方式如下所示：

方式一：

采用循环的方式，删除每一行最后一个数值元素前的非数值元素(非数值元素指isnan函数处理后数值1对应的原矩阵中的元素)，该行剩余元素位置前移，再通过判断矩阵中数值元素的列数，进行矩阵分割，得到节点矩阵和单元矩阵；

方式二：

通过判断原始矩阵的各数值块的列数，识别数值块属于节点矩阵或单元矩阵，再将同一类别的数值块分割拼接，得到节点矩阵和单元矩阵；所述数值块中每行的列数相同且列数的序号一致。采用isnan和sum函数识别数值块类别，获得节点矩阵和单元矩阵的行数和列数，通过矩阵分割成节点矩阵Node和单元矩阵C3D8。

isnan函数处理后为了便于理解，重新赋值，其中1表示节点数据或单元数据等数值数据，0表示非数值数据。然后通过sum函数对各列求和，通过比较各列的数值关系，得到具体的单元数据所在的行列。

将同一类别的数值块进行提取拼接，得到节点矩阵和单元矩阵。

同时，也可以采用其他方式识别各数值块，例如循环的方式；但上述方式，首先通过isnan函数识别非数值元素，再通过sum函数对各列求和，通过比较各列的数值关系，得到具体的单元数据所在的行列，仅进行简单的加减计算，即可得到结果，计算方式简单快速，相较于方式一或采用循环识别数值块的方式，计算效率更高。

步骤S32，扩展节点矩阵，并重构单元矩阵；对单元矩阵内的重复节点重新编号，生成重构单元矩阵，并且扩展节点矩阵包含重新编号的节点；根据原单元矩阵获取表示单元组成面的原面矩阵，根据重构单元矩阵获得表示重构单元组成面的重构面矩阵；

每一单元由6个面、12条线条组成，每一面一般被2个单元共用，每一线条一般被4个单元共用，每一节点一般被8个的单元共用，为了全局嵌入三维Cohesive单元，需将各节点、各线条、各面进行扩展，将单元共用的节点进行区分，使各单元独立。步骤S32包括以下步骤：

步骤S321，扩展节点矩阵，节点矩阵扩展方式为，提取节点矩阵的第一列数据节点编号，将原节点编号扩展为多个不重复的扩展节点编号，使每个单元的节点独立存在，且扩展后的扩展节点对应的三维坐标值为原节点的三维坐标值。

例如，采用矩阵分割和叠加等处理，将节点矩阵的第一列编号矩阵分割出来，并分别加上1000000、2000000、3000000、4000000、5000000、6000000、7000000、8000000、9000000，形成新的编号矩阵Node_1000000、Node_2000000、Node_3000000、Node_4000000、Node_5000000、Node_6000000、Node_7000000、Node_8000000、Node_9000000，将以上的编号矩阵按等列拼接在一起，形成新的编号矩阵Node_No。把原节点矩阵Node中表示节点坐标X轴向、Y轴向、Z轴向坐标的三列(即Node矩阵的第二列、第三列、第四列)提取出来形成坐标矩阵Node_xyz_234，复制形成九个矩阵，并按等列拼接成一个矩阵Node_XYZ，随后将矩阵Node_No和矩阵Node_XYZ按等行拼接在一起，形成行数为原节点矩阵9倍，列数与原节点矩阵一致，编号不同而坐标有重复的新的扩展节点矩阵，记作Inp_3D_Node；本实施例采用节点编号数据递增的方式扩展节点编号，除此之外，还可采用增加不同标志位等其他现有技术扩展节点编号。

步骤S322，重构单元矩阵，单元矩阵的重构方式为，对单元矩阵包含的重复节点进行重新编号，使单元矩阵包含的节点不重复；重复节点的编号方式与节点矩阵的节点编号扩展方式相对应。

例如，本实施例采用以下步骤重构单元矩阵：

步骤S3221，更改单元矩阵C3D8的单元编号列，区分单元编号与节点编号的数据区间；将更改单元编号列后的单元矩阵变形，生成单行或单列的单元矩阵C3D8_1。

例如，切割矩阵，删除单元矩阵Ele_4node的单元编号列，得到第一单元矩阵；并在第一单元矩阵中添加新的单元编号列，新的单元编号列通常添加于第一单元矩阵第一列之前或第一单元矩阵最后一列之后。

步骤S3221中，新的单元编号列所组成的矩阵为过渡矩阵，该过渡矩阵是为方便单元内节点重编号所取，需要跟单元矩阵中的节点编号不重复，一般而言，数值模拟建立的模型中的节点数很少达到上千万，上亿个。因此，为了满足大部分模拟需求，本实施例从100000000开始倒序排序得到新的单元编号列。采用较大值开始的倒序排序的方式是为了使得过渡矩阵中的数字与单元矩阵中的节点编号不重复。除此之外，该过渡矩阵中的新的单元编号列也可取节点编号最大值开始的顺序排列。

优选地，采用reshape函数将更改单元编号列后的单元矩阵变形，生成单行或单列的单元矩阵。

步骤S3222，获取节点编号不重复的节点数目值Unique_C3D8L，建立for循环，根据for循环，查找各节点在单行或单列的单元矩阵中出现的次数及位置，对单行或单列的单元矩阵包含的重复节点进行重新编号，重复节点的编号方式与节点矩阵的节点编号扩展方式相对应；单行或单列的单元矩阵的重复节点重新编号后，通过矩阵分割或变形，转换为与原始单元矩阵C3D8大小一致的重构单元矩阵Inp_3D_C3D8。

优选地，采用unique函数和length函数获取节点编号不重复的节点数目值Unique_C3D8L；采用unique函数去掉已删除单元编号列的单元矩阵中的重复元素，得到单值矩阵Unique_C3D8，采用length函数计算矩阵Unique_C3D8的长度Unique_C3D8L。除此之外，还可根据节点矩阵得到节点数目值Unique_C3D8L，但该方法可能存在误差，根据unique函数和length函数获得单元矩阵实际包含的不重复的节点数目值，数据更准确。所述节点编号不重复的节点数目值Unique_C3D8L用于构建for循环的循环条件，使循环不遗漏节点，并且及时结束循环。

优选地，for循环内采用ismember函数及bwlabel函数协助，以得到重复节点出现次数及位置；ismember函数用于将单行或单列的单元矩阵二值化，以判断单元矩阵中的数据是否为待判断的节点；bwlabel函数用于对二值化后的单元矩阵进行连通域判断，以得到重复节点出现次数及位置。

ismember函数和bwlabel函数通常应用于图像处理，特别是bwlabel函数通常应用于二值图像的连通域判断，通过采用ismember函数及bwlabel函数，简化常规采用for循环以查找各节点在单元矩阵中出现的次数及位置的编程过程，并提高了算法效率。由于每一单元内包含的节点一定不同，且单行或单列的单元矩阵中，相邻单元间通过单元编号隔开，因此通过ismember函数得到的二值化的单元矩阵可以表示某一节点出现的位置，且节点出现的位置一定不相邻，因此可以采用bwlabel函数，通过判断连通域，快速地得到某一位置的节点出现的次数，以得到重复节点出现次数及位置。

S323，原单元矩阵3D_C3D8和重构单元矩阵Inp_3D_C3D8的节点编号列组合，得到原面矩阵C3D80_S6和重构面矩阵C3D81_S6，所述原面矩阵和重构面矩阵用于表示组成六面体的六个面。

具体地，本实施例中将原六面体单元矩阵3D_C3D8的第一列单元编号去除，建立原六面体单元的节点编号矩阵C3D80，同时，去除重新编号后的重构单元矩阵Inp_3D_C3D8的单元编号，新的矩阵命名为C3D81。如图5所示，分别提取矩阵C3D80和矩阵C3D81的第一列、第二列、第三列、第四列建立新的矩阵C3D80_1和矩阵C3D81_1；分别提取矩阵C3D80和矩阵C3D81的第五列、第六列、第七列、第八列建立新的矩阵C3D80_2和矩阵C3D81_2；分别提取矩阵C3D80和矩阵C3D81的第一列、第二列、第六列、第五列建立新的矩阵C3D80_3和矩阵C3D81_3；分别提取矩阵C3D80和矩阵C3D81的第二列、第三列、第七列、第六列建立新的矩阵C3D80_4和矩阵C3D81_4；分别提取矩阵C3D80和矩阵C3D81的第三列、第四列、第八列、第七列建立新的矩阵C3D80_5和矩阵C3D81_5，分别提取矩阵C3D80和矩阵C3D81的第四列、第一列、第五列、第八列建立新的矩阵C3D80_6和矩阵C3D81_6。六个矩阵分别对应单元六个面的节点编号，将以上矩阵分别重新组合C3D80_S6＝[C3D80_1；C3D80_2；C3D80_3；C3D80_4；C3D80_5；C3D80_6]，C3D81_S6＝[C3D81_1；C3D81_2；C3D81_3；C3D81_4；C3D81_5；C3D81_6]，形成新的原面矩阵C3D80_S6和重构面矩阵C3D81_S6。

其中，步骤S4的详细过程如下所示：

S41，根据原面矩阵查找各六面体单元所共有的面，并查找共用的面在原面矩阵和重构面矩阵中的位置，然后查询在原面矩阵和重构面矩阵中用于表示共用的面的节点编号，查询的原画矩阵的节点编号根据工程处理软件的编号规则进行重新排序；查询的重构面矩阵的节点编号用于表示三维Cohesive单元的组成节点，查询的重构面矩阵的节点编号根据原画矩阵中节点的重新排序的规则排序后生成用于添加三维Cohesive单元的三维Cohesive矩阵。

步骤S323所述的重构面矩阵C3D81_S6即为ABAQUS孤立三维模型的所有面的节点编号矩阵，即重构面矩阵C3D81_S6包含构成各三维单元的六个面的节点编号信息。为在单元交界面嵌入Cohesive单元，需判定哪些面为六面体单元的公共面，哪些为非公共面，在公共面添加新的立体单元，以嵌入三维Cohesive单元，再将该边界的所有节点按ABAQUS的INP文件规则重新排列。将原面矩阵中表示六面体单元各面的四列节点编号数字分别相加和相乘，之后再将获得的和与乘积以字符串进行组合，组合后将其重新转化为数字，定义为变量Product_Sum。利用unique函数查找Product_Sum中的单值，建立for循环，通过find函数寻找Product_Sum矩阵中不同变量的位置矩阵Sit_find，输出位置矩阵的大小；若位置矩阵长度为2，则查找在原面矩阵C3D80_S6中Sit_find矩阵第一个数字对应位置处用于表示面的四个节点编号，并存储在循环变量的数字行下的矩阵Cohe001的前四列中。同时，在节点矩阵Inp_3D_Node中查找并输出四个节点的三维坐标位置，新建变量为Cohe4_xyz。以坐标原点为参考点，采用if函数进行判定，使该循环下Cohe001矩阵该行前四列的节点坐标满足闭合和逆时针排列，实现该面四个坐标的重新排序，保存为矩阵Cohe011。查找在原面矩阵C3D80_S6中Sit_find矩阵第二个数字对应位置处用于表示面的四个节点编号，存为矩阵Cohe002。查找在重构面矩阵C3D81_S6中Sit_find矩阵第一个数字对应位置处用于表示面的四个节点编号，存为矩阵Cohe101；查找在重构面矩阵C3D81_S6中Sit_find矩阵第二个数字对应位置处用于表示面的四个节点编号，存为矩阵Cohe102。矩阵Cohe101和矩阵Cohe102根据Cohe011重新排序的规则排序后得到矩阵Cohe111和矩阵Cohe112。具体的，可以根据通过ismember函数获得矩阵Cohe011中的数字在矩阵Cohe001中的位置矩阵Cohe_sit01，按位置矩阵重新排列矩阵Cohe111；通过ismember函数获得矩阵Cohe011中的数字在矩阵Cohe002中的位置矩阵Cohe_sit02，按位置矩阵重新排列矩阵Cohe112。将矩阵Cohe111和矩阵Cohe112拼接成一行，存为COH3D8矩阵，COH3D8矩阵用于非流固耦合单元的全局嵌入；若需要全局嵌入COH3D8P流固耦合模拟的COHESIVE单元，则在矩阵Cohe111内每个数字均加上9000000，存为矩阵Cohe113；将矩阵Cohe111、Cohe112、以及Cohe113拼接在一起，存为COH3D8P矩阵，COH3D8P矩阵用于流固耦合单元的全局嵌入。COH3D8矩阵或COH3D8P矩阵即为本步骤计算获得的三维Cohesive矩阵。

由于节点编号不重复，且立体单元不同面，立体单元不同面的节点编号的和与乘积一般不同，因此上述方法通过计算各面的四各节点编号的和与乘积的方式判断不同单元是否共面，且计算和与乘积与节点编号的顺序无关，进一步降低计算难度。除此之外为了保证计算的精确度，还可以通过判断立体单元不同面是否共同含有3个一致的节点或是否包含2个一致的线段的方式来判断是否不同单元是否共面，该方式计算复杂，耗时，不如计算4个节点编号和与乘积的方法简单效率高。图7至图9给出了计算过程的一个实例，仅做参考。

S42，删除扩展节点矩阵，重构单元矩阵以及三Cohesive矩阵多余的数据部分后，放入模型文件的对应位置替换原节点数据和单元数据，得到新模型文件；

建立单元编号矩阵为从1：1000000的间距为1的列矩阵，并建立长度为1000000，矩阵每一行包含的符号为‘，’。已经处理后的数字矩阵(包括扩展节点矩阵，重构单元矩阵以及Cohesive矩阵)转换为字符串矩阵，通过字符串矩阵的截取和拼接形成新的字符串矩阵INP_3D_Node、INP_C3D8、INP_C3D8P或INP_C3D8；

利用dlmwrite函数将以上三个矩阵分别存入到TXT文本中形成文本：INP_3D_Node.txt、INP_C3D8.txt、INP_C3D8P.txt或INP_C3D8.txt；

将三个文件分别放入到步骤S1所述的inp文件的对应位置替换原来的节点和单元数据形成新的inp文件，另存为：3D_COH.inp；

使用ABAQUS软件打开该inp文件3D_COH.inp，读入孤立网格模型数据，获得全局嵌入三维Cohesive单元后的模型，将所有单元隐藏之后，在mesh模块下使用多余节点删除，即可获得完备的孤立网格模型以建立机械、矿业、航天、土木、隧道、桥梁等专业中的裂纹扩展模型，便于裂纹扩展数值模拟计算。处理后的矩阵包含一些冗余的数据，例如扩展节点矩阵扩展时为了避免扩展数量不足，多扩充的节点数据，这部分数据对应的节点并未出现在重构单元矩阵中，需删除这部分节点数据。

本实施例所述基于MATLAB的三维Cohesive单元全局嵌入方法，根据ABAQUS中inp文件内节点和单元编号的特点，利用EXCEL和MATLAB矩阵实验室，对ABAQUS的inp文件中六面体单元的节点和单元矩阵进行数据处理，建立全局Cohesive单元的节点和编号矩阵，再将节点矩阵和单元矩阵按照ABAQUS的inp文件格式修改为在ABAQUS中能够读取的包含节点和单元数据的inp文件。本发明主要利用MATLAB中的各种函数对ABAQUS中inp文件内的节点和单元矩阵进行数据处理，能够建立全局嵌入三维Cohesive单元的代码库，在短时间将ABAQUS中地质体几何模型实现全局嵌入Cohesive单元，极大提高了ABAQUS的建模速率，且采用外部程序对inp文件的几何数据进行处理，降低了ABAQUS版本变化对全局嵌入Cohesive单元的影响。

本实施例所采用的MATLAB函数仅为本发明具体实施方式的详细说明，可以采用MATLAB中其他功能类似的函数或自编函数；且数据分析软件不限于MATLAB。

以上所述，仅为本发明具体实施方式的详细说明，而非对本发明的限制。相关技术领域的技术人员在不脱离本发明的原则和范围的情况下，做出的各种替换、变型以及改进均应包含在本发明的保护范围之内。

Claims (10)

1.一种三维Cohesive单元全局嵌入方法，其特征在于，包括以下步骤：

S1，采集裂纹信息，通过工程处理软件建立三维的有限元模型，进行网格剖分，并导入实体，处理后的文件保存为模型文件，并输出该模型文件；

S2，通过数据处理软件读取三维的模型文件，提取模型文件中的单元数据及节点数据，并将提取的数据保存为数据文件，并输出该数据文件；

S3，采用数据分析软件读取数据文件记作原始矩阵，原始矩阵通过矩阵分割获得单元矩阵和节点矩阵；分别对单元矩阵和节点矩阵进行处理，将各立体单元共用的节点及面区分开，分别得到扩展节点矩阵和重构单元矩阵；单元矩阵和重构单元矩阵的节点编号列组合，得到原面矩阵和重构面矩阵；

S4，根据原面矩阵获取各六面体单元共用的面，再根据重构面矩阵获得三维Cohesive矩阵；将扩展节点矩阵，重构单元矩阵以及三维Cohesive矩阵放入模型文件的对应位置替换原节点数据和单元数据，得到新模型文件；

S5，工程处理软件读取新模型文件，得到全局嵌入三维Cohesive单元的模型，进行模拟分析。

2.根据权利要求1所述的三维Cohesive单元全局嵌入方法，其特征在于，所述步骤S3包括以下步骤：

步骤S31，采用数据分析软件读取数据文件记作原始矩阵，区分原始矩阵的单元数据和节点数据，获得单元矩阵和节点矩阵；

步骤S32，对节点矩阵进行扩展，并对单元矩阵进行重新构建；对单元矩阵内的重复节点重新编号，生成重构单元矩阵，并且扩展节点矩阵包含重新编号的节点；根据原单元矩阵获取表示单元组成面的原面矩阵，根据重构单元矩阵获得表示重构单元组成面的重构面矩阵。

3.根据权利要求2所述的三维Cohesive单元全局嵌入方法，其特征在于，所述步骤S31采用MATLAB软件作为数据分析软件，首先判定原始矩阵中的非数值元素，排除非数值元素的干扰，通过矩阵分割生成节点矩阵和单元矩阵。

4.根据权利要求3所述的三维Cohesive单元全局嵌入方法，其特征在于，所述步骤S31采用isnan函数和sum函数进行矩阵分割，获得单元矩阵和节点矩阵。

5.根据权利要求2所述的三维Cohesive单元全局嵌入方法，其特征在于，所述步骤S32包括以下步骤：

步骤S321，对节点矩阵进行扩展，节点矩阵的扩展方式为，提取节点矩阵的第一列数据节点编号，将原节点编号扩展为多个不重复的扩展节点编号，使每个单元的节点独立存在，且扩展后的扩展节点对应的平面坐标值为原节点的平面坐标值；

步骤S322，对单元矩阵进行重新构建，单元矩阵的重构方式为，对单元矩阵包含的重复节点进行重新编号，使单元矩阵包含的节点不重复；重复节点的编号方式与节点矩阵的节点编号扩展方式相对应；

S323，单元矩阵和重构单元矩阵的节点编号列组合，得到原面矩阵和重构面矩阵；所述原面矩阵和重构面矩阵用于表示组成六面体的六个面。

6.根据权利要求5所述的三维Cohesive单元全局嵌入方法，其特征在于，所述步骤S322包括以下步骤：

步骤S3221，更改单元矩阵的单元编号列，区分单元编号与节点编号的数据区间；将更改单元编号列后的单元矩阵变形，生成单行或单列的单元矩阵；

步骤S3222，获取节点编号不重复的节点数目值，建立for循环，根据for循环，查找各节点在单行或单列的单元矩阵中出现的次数及位置，对单行或单列的单元矩阵包含的重复节点进行重新编号，重复节点的编号方式与节点矩阵的节点编号扩展方式相对应；单行或单列的单元矩阵的重复节点重新编号后，通过矩阵分割或变形，转换为与原始单元矩阵大小一致的重构单元矩阵。

7.根据权利要求6所述的三维Cohesive单元全局嵌入方法，其特征在于，所述步骤S3222中采用unique函数和length函数获取节点编号不重复的节点数目值Unique_C3D8L；采用unique函数去掉已删除单元编号列的单元矩阵中的重复元素，得到单值矩阵Unique_C3D8，采用length函数计算矩阵Unique_C3D8的长度Unique_C3D8L；所述步骤S3222中for循环内通过ismember函数及bwlabel函数，得到重复节点出现次数及位置；ismember函数用于将单行或单列的单元矩阵二值化，以判断单元矩阵中的数据是否为待判断的节点；bwlabel函数用于对二值化后的单元矩阵进行连通域判断，以得到重复节点出现次数及位置。

8.根据权利要求5所述的三维Cohesive单元全局嵌入方法，其特征在于，所述步骤S4包括以下步骤：

S41，根据原面矩阵查找各六面体单元所共有的面，并查找共用的面在原面矩阵和重构面矩阵中的位置，然后查询在原面矩阵和重构面矩阵中用于表示共用的面的节点编号，查询的原画矩阵的节点编号根据工程处理软件的编号规则进行重新排序；查询的重构面矩阵的节点编号用于表示三维Cohesive单元的组成节点，查询的重构面矩阵的节点编号根据原画矩阵中节点的重新排序的规则排序后生成用于添加三维Cohesive单元的三维Cohesive矩阵；

S42，删除扩展节点矩阵，重构单元矩阵以及三Cohesive矩阵多余的数据部分后，放入模型文件的对应位置替换原节点数据和单元数据，得到新模型文件。

9.根据权利要求8所述的三维Cohesive单元全局嵌入方法，其特征在于，所述步骤S41询共用面的位置的方式为：将原面矩阵中表示六面体单元各面的四列节点编号数字分别相加和相乘，之后再将获得的和与乘积以字符串进行组合，组合后将其重新转化为数字，定义为变量Product_Sum；利用unique函数查找Product_Sum中的单值，建立for循环，通过find函数寻找Product_Sum矩阵中不同变量的位置矩阵Sit_find，若位置矩阵长度为2，则其对应的单元共面，共用面为位置矩阵指示的面。

10.根据权利要求8所述的三维Cohesive单元全局嵌入方法，其特征在于，所述三维Cohesive单元为流固耦合单元或非流固耦合单元。

Images