二次开发黏聚裂纹子程序

2023-06-27 理论教育版权反馈

【摘要】：有限元法是目前发展比较完善的数值分析方法，目前常见的通用软件ABAQUS、LS－DYNA、ANSYS 提供了很多应用程序接口和用户子程序接口，方便用户根据专业问题定义自己所需的模型，并通过二次开发建立所需要的本构模型。黏聚裂纹子程序是通过LS－DYNA 二次开发实现的。根据第3 章介绍的黏聚裂纹模型理论，有限元程序在每一时间步提供的应变增量为Δεa，在每次计算前需要提取前一步存储的历史变量。表4－1UMAT 中的子程序

对PBX 这类准脆性材料的断裂过程进行分析，可以发现其失效是微裂纹逐渐会聚成为宏观裂纹并断裂的过程。黏聚裂纹模型［1，2］所需材料参数简单，物理意义明确，可以描述材料任意位置拉伸裂纹的起裂或扩展，以及裂纹扩展导致的材料软化行为。嵌入式黏聚裂纹模型基于其基本思想，通过单元形函数的扩展，将黏聚裂纹模型直接作用于计算网格内部，避免了处理裂纹尖端奇异的困难。

有限元法是目前发展比较完善的数值分析方法，目前常见的通用软件ABAQUS、LS－DYNA、ANSYS 提供了很多应用程序接口和用户子程序接口，方便用户根据专业问题定义自己所需的模型，并通过二次开发建立所需要的本构模型。用户材料子程序（User defined material，UMAT）是提供的自定义材料程序接口，用户可以通过此接口扩展有限元计算的适用范围（图4－1）。黏聚裂纹子程序是通过LS－DYNA 二次开发实现的。

有限元计算（增量方法）的基本问题是：已知第n 步的结果（应力、应变等）σn、εn，然后给出一个应变增量dεn＋1，计算σn＋1，UMAT 的主要任务就是完成这一计算。根据第3 章介绍的黏聚裂纹模型理论，有限元程序在每一时间步提供的应变增量为Δεa，在每次计算前需要提取前一步存储的历史变量。如图4－1 所示，为在主程序中定义的crack_flag （判定单元中裂纹是否存在）、 pagenumber_ebook=76,pagenumber_book=56 （裂纹历史最大张开位移）、b＋（形函数梯度矢量）、n （裂纹法向矢量）、w （裂纹张开位移矢量）。这些变量在每一步结束后存储于历史变量中，以便于下一步计算前提取。

pagenumber_ebook=77,pagenumber_book=57

图4－1　子程序流程图

在初始阶段，单元内部没有裂纹出现，则crack_flag＝0，形函数梯度矢量b＋和裂纹法向矢量n 未定义，单元变形为连续体变形，所以应力σ 可以直接通过连续体材料本构进行计算。当最大主应力σI超过拉伸强度ft时，黏聚裂缝将在垂直于最大主应力的方向生成，此时crack_flag＝1。这时，嵌入的裂纹将连续体分隔开来，Δεa分为两部分：连续体单元的应变增量Δεc和裂纹张开位移引起的变形增量。为了计算裂纹张开位移，根据n （与最大主应力特征矢量方向相同），开始计算b＋。当裂纹嵌入后，裂纹法向矢量n 将保持恒定，在每一个时间步将计算裂纹张开位移矢量w，在这里需要提取前一步的w，计算裂纹张开位移的增量Δw，随后计算Δεc，从而可以计算得到应力增量。考虑到载荷卸载条件，对比裂纹历史最大张开位移 pagenumber_ebook=77,pagenumber_book=57 和计算得到的裂纹张开位移的模｜w｜，如果｜w｜＞，则将｜w｜储存为的值。

为了方便编译和读取，UMAT 中划分为多个子程序，功能列入表4－1 中。其中，Sub.1 中计算连续体单元的本构方程主要有线弹性本构和弹塑性本构；Sub.5 中有多种软化曲线的模式，包括矩形、线性、指数型软化曲线。由于，形函数梯度及裂纹方向是根据各单元对应节点的坐标得来的，所以在Sub.8中提供单元的节点坐标信息。为了获取当前正在计算的单元编号，在UMAT中增加了新的参量，并相应地在调用位置进行声明。声明传递参量后，在每个循环中，由内部函数nelmntid 读取当前的单元号。根据单元号读取外部文件中，获取对应的节点坐标信息，并将其存入历史变量中，用于后续计算。在UMAT 中需要用户指定的材料参数，如体积模量、剪切模量、拉伸强度、断裂能、软化曲线及材料本构选择等。

表4－1　UMAT 中的子程序

pagenumber_ebook=78,pagenumber_book=58

二次开发黏聚裂纹子程序

相关推荐