(12)发明专利申请
(10)申请公布号 CN 112749501 A(43)申请公布日 2021.05.04
(21)申请号 202110036380.4(22)申请日 2021.01.12
(71)申请人 哈尔滨工业大学
地址 150001 黑龙江省哈尔滨市南岗区西
大直街92号(72)发明人 果立成 沈日麟 郝留磊 于红军
申振 (74)专利代理机构 哈尔滨龙科专利代理有限公
司 23206
代理人 高媛(51)Int.Cl.
G06F 30/23(2020.01)G06F 30/17(2020.01)G06F 111/10(2020.01)G06F 119/14(2020.01)
权利要求书2页 说明书5页 附图6页
(54)发明名称
一种基于幂指数型准则的断裂行为仿真相场法
(57)摘要
本发明公开了一种基于幂指数型准则的断裂行为仿真相场法，所述方法如下：一、引入Mα和Mβ量化不同的能量对裂纹演化的影响，基于变分原理推导相场法的控制方程以及有限元离散格式；二、通过Fortran语言编写相场单元程序、计算执行文件和输出命令文件，并定义边界条件、载荷信息和材料属性信息；三、在ABAQUS中建立有限元模型并划分网格，导出相应的单元和节点信息，生成模型信息文件；四、利用FEAP快速求解器求解，输出载荷、位移、应力、应变和相场等信息，并将求解结果导入ParaView中得到可视化结果。本发明的方法不仅过程简单，而且具有良好的适应性与稳定性，可开发成商用有限元程序以灵活适应问题的变化。
CN 112749501 ACN 112749501 A
权 利 要 求 书
1/2页
1.一种基于幂指数型准则的断裂行为仿真相场法，其特征在于所述断裂相场模拟方法包括如下步骤：
步骤一：基于相场理论，从裂纹面拓扑表征的角度出发，引入幂指数参数Mα和Mβ来量化不同的能量对裂纹演化的影响，并基于变分原理推导相场法的控制方程以及相应的有限元离散格式；
步骤二：通过Fortran语言编写基于幂指数模型的相场单元程序、计算执行文件和输出命令文件，并在计算执行文件中定义边界条件、载荷信息和材料属性信息；
步骤三：在ABAQUS中建立有限元模型并划分网格，在预期的裂纹扩展区域进行细化，并导出相应的单元和节点信息，生成模型信息文件；
步骤四：利用FEAP快速求解器实现对I型和II型混合型断裂问题的求解，采用牛顿插值的方法求解非线性方程，当计算满足收敛准则时进入下一个时间步，输出载荷、位移、应力、应变和相场分布信息，并将求解结果导入至ParaView后处理软件中得到可视化结果。
2.根据权利要求1所述的基于幂指数型准则的断裂行为仿真相场法，其特征在于所述步骤一中，基于变分原理对相场法的控制方程进行推导的推导结果如下：
其中：为哈密顿微分算子，Δ为拉普拉斯算子，mc为退化函数，σ为主应力张量，n为主应力方向张量，b为系统体力，t为作用在系统表面上的面力，l0为长度尺度，u为位移场，为边界上位移场的值，
为整个区域的边界，
为力的边界，
为位移场的边界，c为
相场，Δ为拉普拉斯算子，为I型断裂历史场，释放率，
为II型临界能量释放率。
提出球量‑偏量分割法和主应变方向空间分割法区分
为II型断裂历史场，为I型临界能量
3.根据权利要求1所述的基于幂指数型准则的断裂行为仿真相场法，其特征在于所述步骤一中，针对历史变量
拉伸和压缩状态。
4.根据权利要求3所述的基于幂指数型准则的断裂行为仿真相场法，其特征在于对于球量‑偏量分割法，历史变量
和
的定义如下：
其中，
为球量部分应变能密度，εdev为偏应变，Θ为体积应变，<Θ>+的定义为
2
CN 112749501 A
权 利 要 求 书
μ为剪切模量，t为分析步，S表示位于[1,t]区间内的分析步。
2/2页
5.根据权利要求3所述的基于幂指数型准则的断裂行为仿真相场法，其特征在于对于主应变方向空间分割法，历史变量
和
的定义如下：
其中，λ为拉梅常数，μ为剪切模量，ε为主应变张量。
6.根据权利要求1所述的基于幂指数型准则的断裂行为仿真相场法，其特征在于所述步骤一中，有限元离散格式为：
其中，Ru为位移场的余量，Rc为相场的余量，Nu为位移场的标准双线性形函数，Nc为相场的标准双线性形函数，Bu为位移场形函数的梯度，Bc为相场形函数的梯度。
3
CN 112749501 A
说 明 书
一种基于幂指数型准则的断裂行为仿真相场法
1/5页
技术领域
[0001]本发明属于计算机模拟技术领域，涉及一种断裂相场模拟方法，具体涉及一种基于幂指数型准则的断裂行为仿真相场法。
背景技术
[0002]在宏观尺度下，有限元模拟是一种非常重要的模拟方法，其对指导工业加工、工程建设等具备重要的指导意义。其中，预测材料的损伤、破坏、断裂行为是有限元模拟的一个非常关键的应用。以往对断裂行为的模拟多采用离散的方法如扩展有限元(XFEM)，然而这类方法需要数值追踪物理上不连续的裂纹表面，无法处理多裂纹复杂相互作用问题，如裂纹融合等。为了克服这些困难，断裂相场法被提出。相场法引入了一个标量场d∈[0,1]，其中0代表完好的材料状态，1代表完全的损伤状态即裂纹。基于此方法，相场法基于全局能量最小化原理获得全新条件下的裂纹面拓扑，避免了追踪间断的裂纹表面，可以有效分析复杂的多裂纹相互作用问题，如分叉、汇合和起裂扩展等。[0003]相场法由变分形式的Griffith理论推导得到，Griffith临界能量释放率是相场法中的重要参数。然而在传统相场法中，仅考虑了I型临界能量释放率，在两种情况下可以取得较好的分析结果：[0004](1)I型和II型断裂能接近的材料的断裂分析；[0005](2)I型和II型断裂能相差较大但处于I型裂纹扩展主导的情况。[0006]当材料I型和II型断裂能相差较大，且失效并不是I型的裂纹主导时，就需要综合考虑I型和II型断裂能对断裂行为的影响。针对材料的裂纹开裂方式可以分为三种：I型面内张开型裂纹、II型面内剪切型裂纹和III型面外剪切型裂纹。在日常服役环境下，材料往往受到多种类型载荷的共同作用，呈现出混合型裂纹扩展的形态。因此，开发一种可同时考虑I型和II型混合型断裂的新型相场法对于相场方法的发展具有十分重要的意义。发明内容
[0007]本发明的目的是提供一种基于幂指数型准则的断裂行为仿真相场法，通过引入I型和II型临界能量释放率并考虑不同类型临界能量释放率的贡献来改进传统相场法，其本质是两项裂纹驱动力
和
分别引入不同的比例因子，通过改变两项裂纹驱动力的相
对比值，最终导致裂纹扩展方向的改变，从而实现对混合型裂纹扩展问题的模拟。[0008]本发明的目的是通过以下技术方案实现的：
[0009]一种基于幂指数型准则的断裂行为仿真相场法，包括如下步骤：[0010]步骤一：基于相场理论，从裂纹面拓扑表征的角度出发，引入幂指数参数Mα和Mβ来量化不同的能量对裂纹演化的影响，并基于变分原理推导相场法的控制方程以及相应的有限元离散格式；[0011]步骤二：通过Fortran语言编写基于幂指数模型的相场单元程序、计算执行文件和
4
CN 112749501 A
说 明 书
2/5页
输出命令文件，并在计算执行文件中定义边界条件、载荷信息和材料属性信息；[0012]步骤三：在ABAQUS中建立有限元模型并划分网格，在预期的裂纹扩展区域进行细化，并导出相应的单元和节点信息，生成模型信息文件；[0013]步骤四：利用FEAP快速求解器实现对I型和II型混合型断裂问题的求解，采用牛顿插值的方法求解非线性方程，当计算满足收敛准则时进入下一个时间步，输出载荷、位移、应力、应变和相场分布等信息，并将求解结果导入至ParaView后处理软件中得到可视化结果。[0014]相比于现有技术，本发明具有如下优点：
[0015]
1、针对两项裂纹驱动力和
本发明分别引入不同的幂指数Mα和Mβ来综合考
虑I型和II型临界能量释放率的权重对裂纹扩展速度和方向的影响，从而实现对混合型裂纹扩展问题的模拟，弥补了传统相场法仅能分析材料处于I型裂纹扩展主导的不足。[0016]2、通过改变模型信息、边界条件、载荷信息、材料属性信息、能量分割方式、幂指数参数Mα和Mβ，可以实现对不同材料在不同受力状态下的I型和II型混合型断裂问题的模拟。[0017]3、本发明的方法不仅过程简单，而且具有良好的适应性与稳定性，可开发成商用有限元程序以灵活适应问题的变化。
附图说明
[0018]图1为基于幂指数型准则的断裂行为仿真相场法的流程框图；[0019]图2为含斜裂纹受压岩石试样的几何尺寸、边界条件、载荷和网格设置示意图；[0020]图3为采用幂指数模型相场法模拟的含裂纹受压岩石试样的裂纹演化历程图；[0021]图4为采用经典相场法模拟的含裂纹受压岩石试样的裂纹形貌图；[0022]图5为含裂纹受压岩石试样在试验中的典型裂纹形貌示意图；[0023]图6为含倾斜裂纹的对称半圆形巴西圆盘试样的几何尺寸、边界条件、载荷以及网格设置示意图；
[0024]图7为采用幂指数模型相场法模拟的含倾斜裂纹的对称半圆形巴西圆盘试样的数值裂纹路径与试验结果对比图，(a)试验结果
(d)模拟结果(h)试验结果(l)模拟结果
具体实施方式
[0025]下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。
[0026]本发明提供了一种基于幂指数型准则的断裂行为仿真相场法，如图1所示，包括如下步骤：
[0027]步骤一：基于相场理论，从裂纹面拓扑表征的角度出发，引入两个新的参数Mα和Mβ
5
(b)试验结果(f)模拟结果(j)模拟结果
(c)试验结果(g)试验结果(k)模拟结果
(e)模拟结果(i)试验结果
CN 112749501 A
说 明 书
3/5页
来量化不同的能量对裂纹演化的影响。幂指数Mα和Mβ是表征材料内部各种缺陷的细观材料参数，这些细观缺陷包括微裂纹、微孔洞、位错和滑移等，幂指数Mα和Mβ可以综合考虑I型和II型临界能量释放率的权重对裂纹扩展速度和方向的影响，在实际应用中需要对Mα和Mβ进行标定。
[0028]基于变分原理对相场法的控制方程进行推导，推导结果如下：
[0029][0030][0031][0032][0033][0034]
其中：为哈密顿微分算子，Δ为拉普拉斯算子，mc为退化函数，σ为主应力张量，n
为整个区域的边界，
为力的边界，
为位移场的边
为主应力方向张量，b为系统体力，t为作用在系统表面上的面力，l0为长度尺度，u为位移场，为边界上位移场的值，
界，c为相场，Δ为拉普拉斯算子，为I型断裂历史场，为II型断裂历史场，为I型临界能量释放率，为II型临界能量释放率。
[0035]
为了避免材料在压缩状态下产生损伤，以及材料断裂后因裂纹面闭合产生的不合
提出了两类能量分割方式来解决这一问题，
理的材料侵入现象，针对历史变量
即球量‑偏量分割法和主应变方向空间分割法。
[0036][0037][0038][0039][0040][0041][0042]
对于球量‑偏量分割法，历史变量和的定义如下：
对于主应变方向空间分割法，历史变量和的定义如下：
其中，
dev
为球量部分应变能密度，ε为偏应变，Θ为体积应变，<Θ>+的定义为
λ为拉梅常数，μ为剪切模量，ε为主应变张量，t为分析步，S表示位于[1,t]区
间内的分析步。
[0043]在考虑拉伸压缩异性的基础上，基于幂指数模型相场法的有限元离散格式为：
[0044]
6
CN 112749501 A
说 明 书
4/5页
[0045]
其中，Ru为位移场的余量，Rc为相场的余量，Nu为位移场的标准双线性形函数，Nc为
相场的标准双线性形函数，Bu为位移场形函数的梯度，Bc为相场形函数的梯度。[0047]步骤二：通过Fortran语言编写基于幂指数模型的相场单元程序、计算执行文件和输出命令文件，相场单元程序命名为elmt15，计算执行文件命名为IZ1L4，输出文件命名为solve.Z1L4，并在计算执行文件中定义边界条件、载荷信息和材料属性信息。[0048]步骤三：在ABAQUS中建立有限元模型并划分网格，为了能够准确地捕捉裂纹形貌，在预期的裂纹扩展区域进行了细化，并导出相应的单元和节点信息，生成模型信息文件，文件命名为ISquarel。[0049]步骤四：将相场单元程序elmt15、计算执行文件IZ1L4、输出命令文件solve.Z1L4和模型信息文件ISquarel导入至装有FEAP的计算机中，然后输入相应的命令进行求解计算。利用FEAP快速求解器实现对I型和II型混合型断裂问题的求解，采用牛顿插值的方法求解非线性方程，当计算满足收敛准则时进入下一个时间步，输出载荷、位移、应力、应变和相
相比于ANSYS等商业有限元软件，开源FEAP有限元提供了较好的二次开发平场分布等信息。
台，并将求解结果导入至ParaView后处理软件中得到可视化结果。[0050]通过改变模型信息、边界条件、载荷信息、材料属性信息、能量分割方式、幂指数参数Mα和Mβ，可以实现对不同材料在不同受力状态下的I型和II型混合型断裂问题的模拟。[0051]为了说明本发明上述方案的性能，下面结合两个算例来进一步描述。[0052]算例一：含斜裂纹受压岩石试样
[0053]含斜裂纹受压岩石试样的几何尺寸、边界条件、载荷和网格设置如图2所示，模型的边长为L＝102mm，在模型中心处预制了长和宽分别为20.32mm和0.204mm。初始裂纹与全局x轴的夹角为45°。当前模型采用位移载荷控制，模型上边界施加了向下的位移载荷u*，位移增量为Δu＝2.0e‑4mm。模型下边界沿y方向固定，且模型左下角度沿x方向固定以避免在准静态加载时由于刚体运动导致的矩阵奇异。为了能够正确地捕捉裂纹形貌，网格在预期的裂纹扩展区域进行了细化。细化区域的有效网格尺寸大约为0.1mm，相场法的长度尺度l0选为三个单元的长度，即0.3mm。模型大约包含了188685个完全积分的四边形单元，共有188776个节点，即566328个自由度。材料的杨氏模量E＝36.2GPa，泊松比υ为0.21。I型断裂能
而II型断裂能为
即110J/m2。设定Mα＝Mβ＝1。采用主
[0046]
应变空间分解法来计算模型的裂纹演化过程。
[0054]图3显示了采用幂指数模型相场法模拟的含裂纹受压岩石试样的裂纹演化过程，
由于拉伸主应变的存在，模型依然会出现裂纹的可以看到，尽管模型处于压缩载荷作用下，
扩展。具体的裂纹扩展过程为：首先出现反对称翼裂纹的扩展演化，扩展一定距离后，翼裂纹停止扩展，而反对称的次生裂纹开始萌生扩展，并导致模型的最终失效。这与文献中实验结果所显示的断裂机制一致，如图5所示，证明了当前模型可以有效考虑岩石混合型裂纹的扩展。为了进一步说明本发明方法的有效性，图4显示了采用经典相场法获得的裂纹形貌图，可以看到，经典相场法模型所预测的裂纹形貌与实验结果有明显不同。[0055]算例二：含倾斜裂纹的对称半圆形巴西圆盘试样
7
CN 112749501 A[0056]
说 明 书
5/5页
含倾斜裂纹的对称半圆形巴西圆盘试样的几何尺寸、边界条件、载荷以及网格设
置如图6所示，模型半径为R＝50mm，对称线上预制了长度为l的初始裂纹。模型对称线顶端的位置施加了向下的位移载荷u*。模型下端由两个对称点支撑，距离对称线的距离分别为S。在当前研究中，取S/R＝0.43和l/R＝0.3保持固定不变。当
的时候，模型为纯II型
开裂，因此初始裂纹和对称线的夹角的分别设置为0°,10°,20°,30°,40°,50°，可以在给定载荷条件下获得不同混合比初始裂纹。为了能够正确地捕捉裂纹形貌，网格在预期的裂纹
相场法的长度尺度l0选择三扩展区域进行了细化。细化区域的有效网格尺寸大约为0.1mm，
个单元的长度，即0.3mm。对于不同初始偏转角，模型细化区域有所差别，但各个模型均包含了4万以上个结构化的四边形单元。此处通过几何划分采取结构化网格，从而保证所有的细化网格尺寸基本一致。材料的杨氏模量E＝2790MPa，泊松比υ＝0.35。在这里设定Mα＝Mβ＝1，
[0057]
图7显示了具有不同初始裂纹偏转角的模型的数值结果和实验结果。从图中可以
在不断增大，而裂纹偏转角
看到，随着初始裂纹偏转角的增大，裂纹尖端的混合比
也在不断增大。对于固定的初始裂纹偏转角，数值结果和实验结果完美吻合。这一算例表明，基于幂指数型准则的断裂行为仿真相场法可以准确描述材料混合型裂纹扩展问题。[0058]以上两个算例验证了本发明上述方案的正确性以及高效性。
8
CN 112749501 A
说 明 书 附 图
1/6页
图1
9
CN 112749501 A
说 明 书 附 图
2/6页
图2
10
CN 112749501 A
说 明 书 附 图
3/6页
图3
11
CN 112749501 A
说 明 书 附 图
4/6页
图4
图5
12
CN 112749501 A
说 明 书 附 图
5/6页
图6
13
CN 112749501 A
说 明 书 附 图
6/6页
图7
14
因篇幅问题不能全部显示，请点此查看更多更全内容