本发明属于道路工程数据处理技术领域，尤其涉及一种基于vcct确定不同裂纹长度的能量释放率的优化方法。

背景技术：

虚拟裂纹扩展分析技术(vcct)是基于裂纹扩展时释放的能量等于将裂纹闭合一个扩展增量所需的功的假设的一种断裂力学分析方法，该方法假设虚拟裂纹尖端后面的展开位移和实际裂纹尖端后面的张开位移近似相等，通过虚拟裂纹扩展分析技术可获得应力或应变的精确解，但基于虚拟裂纹扩展分析技术现有的abaqus中的一个模型只能获得单一裂纹长度对应的能量释放率，无法在单个模型中获得不用裂纹长度对应的能量释放率，若要获得不同裂纹长度对应的能量释放率只能对不同裂纹长度进行建模，如计算n中裂纹长度对应的能量释放率，需要建立n个模型，这样大大增大了建模工作量和计算量，且对计算机的配置要求高。

因此，亟需一种基于虚拟裂纹扩展分析技术确定不同裂纹长度的能量释放率的优化方法。

技术实现要素：

有鉴于此，本发明提供一种基于vcct确定不同裂纹长度的能量释放率的优化方法以解决现有技术的不足。

本发明提供一种基于vcct确定不同裂纹长度的能量释放率的优化方法，其特征在于：包括以下步骤：

s1：在abaqus中建立复合式路面结构模型；在复合式路面模型上预设裂纹和裂纹参数，然后将所述复合式路面模型进行网格划分，其中，所述裂纹参数包括裂纹扩展方向、时间增量步和牵引开裂值；

s2：从abaqus中提取所述裂纹扩展面的粘结节点并将所述粘结节点按照裂纹扩展路径重新排列，形成裂纹扩展节点矩阵；

s3：根据所述裂纹扩展节点矩阵，在abaqus中计算不同裂纹长度对应的能量释放率。

进一步，所述步骤s2包括以下步骤：

s21：将位于预设裂纹两侧的粘结节点分别放置于l矩阵和r矩阵；所述粘结节点至少包括节点坐标和节点编号；

s22：将所述l矩阵按照列进行升序排列形成新的矩阵l1；将所述r矩阵按照列进行升序排列形成新的矩阵r1；

s23：将矩阵l1中的第三列赋值给矩阵nl，将矩阵r1的第三列赋值给矩阵nr；即将l1矩阵的节点编号幅值给矩阵nl，将r1矩阵的节点编号幅值给矩阵nr，即矩阵nl和矩阵nr为节点编号；

s24：定义列向量n1和n2:

s241：判断w与len的大小，其中，w为中间变量，其初始值为1，len表示矩阵nl中的元素个数或nr中的元素个数，若w小于等于len，进一步判断w除row的余数是否是1，其中，row＝len/col，len表示矩阵nl中的元素个数，col表示节点列数，若余数不是1，进入步骤s242，若余数是1，则进入步骤s243,若w大于len,则退出；

s242：n1[k]＝nl[w],n2[k]＝nr[w]，其中，k为列向量n1和n2的行数，w表示阵列nl和阵列nr的行数；

s243：w＝w+1,k＝k+1；

s244：重复步骤s241至步骤s243，直至w大于len退出；

s25：定义列向量n3和n4；

s251：判断j与len的大小，其中，j为中间变量其初始值为1，len表示矩阵nl中的元素个数或nr中的元素个数，若j小于等于len，进一步判断j除row的余数是否是0，其中，row＝len/col，len表示矩阵nl中的元素个数或nr中的元素个数，col表示节点列数，若余数不是0，进入步骤s252，若余数是0，则进入步骤s253,若j大于len,则退出；

s252：n3[n]＝nl[j],n2[n]＝nr[j]，其中，n为列向量n3和n4的行数，j表示阵列nl和阵列nr的行数；

s253：j＝j+1,n＝n+1；

s254：重复步骤s251至步骤s253，直至w大于len退出；

s26：输出裂纹扩展节点矩阵node,node＝[n1n2n3n4]；

s27：将裂纹扩展节点矩阵node作为自定义单元写入aqaqus的inp文件中，形成修改后的inp文件。

进一步，所述裂纹扩展节点矩阵的一行的四个节点表示一个哑结点断裂单元的四个节点。

进一步，其中所述l矩阵和r矩阵均为多行三列矩阵，从左至右依次为，第一列为节点的x轴坐标，第二列为节点的y轴坐标，第三列为节点的编号。

进一步，所述能量释放率包括ⅰ型能量释放率和ⅱ型能量释放率。

进一步，所述不同裂纹长度对应的能量释放率在abaqus中采用如下方法计算：

s31：从inp文件读入自定义单元的第i行的节点，即裂纹扩展节点矩阵node的第i行的节点，其中,i为中间变量，其初始值为1，形成弹性刚度矩阵amatrx，并更新右手端矢量rhs；

s32：计算哑结点断裂单元的能量释放率，并赋值给状态变量矩阵svars，其中，所述状态变量矩阵svars包括svars(1)和svars(2)；

s33：判断分析步时间是否小于等于分析步一，即step-1的结束时间，若是，则弹性刚度矩阵amatrx元素数值不变，裂纹不发生扩展，进入步骤s34；若否，进一步判断状态变量svars的元素svars(2)或svars(1)是否大于预设的牵引开裂临界值，若是，则弹性刚度矩阵amatrx元素全部赋值为0，裂纹扩展，裂纹长度增长△a，进入步骤s34；若否，弹性刚度矩阵amatrx元素数值不变，裂纹不发生扩展，进入步骤s34；

s34：判断当前分析步时间是否大于step-2结束时间，其中，step-2结束时间表示裂纹扩展的计算总时间，若否，则分析步时间增加一个时间增量步，进入步骤s31；若是，退出。

进一步，所述步骤s33还包括将所述状态变量矩阵svars幅值给abaqus状态变量sdv1和sdv2，并输出sdv1和sdv2。

本发明的有益效果：本发明将裂纹扩展面粘结节点按预设裂纹扩展路径进行排序，以及裂纹扩展的条件判断和每个裂纹扩展增量的能量释放率的计算，实现在abaqus的单个模型中对不同裂纹长度对应的能量释放率的计算，具有建模工作量小、计算量小和精度高的有益技术效果。

附图说明

下面结合附图和实施例对本发明作进一步描述：

图1为生成自定义单元matlab程序流程图。

图2为计算不同裂纹长度对应的能量释放率的流程图。

具体实施方式

以下结合说明书附图对本发明做出进一步的说明：

本发明提供一种基于vcct确定不同裂纹长度的能量释放率的优化方法，其特征在于：包括以下步骤：

s1：在abaqus中建立复合式路面结构模型；在复合式路面模型上预设裂纹和裂纹参数，然后将所述复合式路面模型进行网格划分，其中，所述裂纹参数包括裂纹扩展方向、时间增量步和牵引开裂值；

s2：从abaqus中提取所述裂纹扩展面的粘结节点并将所述粘结节点按照裂纹扩展路径重新排列，形成裂纹扩展节点矩阵；

s3：根据所述裂纹扩展节点矩阵，在abaqus中计算不同裂纹长度对应的能量释放率。

其中，vcct表示virtualcrackclosuretechnique的简称，即虚拟裂纹扩展分析技术或虚拟裂纹闭合法的简称。

通过上述技术方案，对裂纹扩展面粘结节点按预设裂纹扩展路径进行排序，以及裂纹扩展的条件判断和每个裂纹扩展增量的能量释放率的计算，实现在abaqus的单个模型中对不同裂纹长度对应的能量释放率的计算，具有建模工作量小、计算量小和精度高的有益技术效果。解决了现有abaqus中，计算不同裂纹长度对用个的能量释放率需要重复建立多个模型，以及计算量大的问题，实现了在单个模型中计算不同裂纹长度对应的能量释放率，即实现了在abaqus中模拟预设裂纹的连续扩展。

在本实施例中，所述步骤s2包括以下步骤：

s21：将位于预设裂纹两侧的粘结节点分别放置于l矩阵和r矩阵；所述粘结节点至少包括节点坐标和节点编号；所述节点坐标可以是平面二维坐标也可以是立体的三维坐标，所述节点编号是abaqus程序对模型内的节点的编号；其中所述l矩阵和r矩阵均为多行三列矩阵，其中所述l矩阵和r矩阵均为多行三列矩阵，从左至右依次为，第一列为节点的x轴坐标，第二列为节点的y轴坐标，第三列为节点的编号。

s22：将所述l矩阵按照列进行升序排列形成新的矩阵l1；将所述r矩阵按照列进行升序排列形成新的矩阵r1；

s23：将矩阵l1中的第三列赋值给矩阵nl，将矩阵r1的第三列赋值给矩阵nr；

s24：定义列向量n1和n2:

s241：判断w与len的大小，其中，w为中间变量，其初始值为1，len表示矩阵nl中的元素个数或nr中的元素个数，若w小于等于len，进一步判断w除row的余数是否是1，其中，row＝len/col，len表示矩阵nl中的元素个数，col表示节点列数，若余数是1，进入步骤s242，若余数不是1，则进入步骤s243,若w大于len,则退出；

s242：n1[k]＝nl[w],n2[k]＝nr[w]，其中，k为列向量n1和n2的行数，w表示阵列nl和阵列nr的行数；

s243：w＝w+1,k＝k+1；

s244：重复步骤s241至步骤s243，直至w大于len退出；

s25：定义列向量n3和n4；

s251：判断j与len的大小，其中，j为中间变量其初始值为1，len表示矩阵nl中的元素个数或nr中的元素个数，若j小于等于len，进一步判断j除row的余数是否是0，其中，row＝len/col，len表示矩阵nl中的元素个数或nr中的元素个数，col表示节点列数，若余数是0，进入步骤s252，若余数不是0，则进入步骤s253,若j大于len,则退出；s252：n3[n]＝nl[j],n2[n]＝nr[j]，其中，n为列向量n3和n4的行数，j表示阵列nl和阵列nr的行数；

s253：j＝j+1,n＝n+1；

s254：重复步骤s251至步骤s253，直至w大于len退出；

s26：输出裂纹扩展节点矩阵node,node＝[n1n2n3n4]；

s27：将裂纹扩展节点矩阵node作为自定义单元写入aqaqus的inp文件中，形成修改后的inp文件。

其中，本步骤的计算本领域技术人员可根据实际选择实现方式，如手动计算、或通过计算机实现，在本实施例中，通过matlab软件来实现。

通过上述技术方案，将裂纹面的粘结节点从复合式道路模型中读取出来，并根据预设的裂纹扩展路径重新排序，使重新排序后矩阵的符合裂纹尖端的哑结点断裂单元的节点排序。

在本实施例中，所述裂纹扩展节点矩阵的一行的四个节点表示一个哑结点断裂单元的四个节点。通过上述技术方案可为后续裂纹的连续扩展提供哑结点断裂单元的扩展路径。

在本实施例中，所述能量释放率包括ⅰ型能量释放率和ⅱ型能量释放率。所述ⅰ型能量释放率表示垂直方向的拉伸裂纹驱动力，所述ⅱ型能量释放率表示水平方向的剪切裂纹驱动力。

在本实施例中，所述不同裂纹长度对应的能量释放率在abaqus中采用如下方法计算：

s31：从inp文件读入自定义单元的第i行的节点，即裂纹扩展节点矩阵node的第i行的节点，其中,i为中间变量，其初始值为1，形成弹性刚度矩阵amatrx，并更新右手端矢量rhs；

s32：计算哑结点断裂单元的能量释放率，并赋值给状态变量矩阵svars，其中，所述状态变量矩阵svars包括svars(1)和svars(2)；

所述ⅰ型能量释放率gⅰ和ⅱ型能量释放率gⅱ采用如下方法计算：

其中，fx、fy分别为节点x轴和y轴方向的节点力，分别为相邻节点编号x轴和y轴方向的位移差，δa为裂纹扩展长度，b表示裂纹厚度；

s33：判断时间增量步是否小于等于分析步一，即step-1的结束时间，此步骤即是判断负载是否加载完成，若是，即负载尚未加载完成，则弹性刚度矩阵amatrx元素数值不变，裂纹不扩展，进入步骤s34；若否，进一步判断状态变量svars的元素svars(2)或svars(1)是否大于预设的牵引开裂临界值，若是，则弹性刚度矩阵amatrx元素全部赋值为0，即弹性刚度矩阵之间以没有牵引力，裂纹扩展，裂纹长度增长△a，进入步骤s34；若否，弹性刚度矩阵amatrx元素数值不变，裂纹不发生扩展，进入步骤s34；

s34：判断当前分析步时间是否大于step-2结束时间，其中，step-2结束时间表示裂纹扩展的计算总时间，若否，则分析步时间增加一个时间增量步，进入步骤s31；若是，退出。

进一步，所述步骤s33还包括将所述状态变量矩阵svars幅值给abaqus状态变量sdv1和sdv2，并输出sdv1和sdv2。通过上述技术方案，将机器可读但用户无法获得的能量释放率转换为机器和用户均可识别的能量释放率。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。