《码上仿真:Abaqus Python API》Vol.14 ContourIntegral 在裂纹尖端画“禁区”

在经典力学里，应力等于力除以面积。但当你面对一个无限尖锐的裂纹时，裂尖的面积趋近于 0，这意味着应力会趋向于无穷大。

这就是数学上的"奇异性（Singularity）"，它就像物理学里的黑洞，让常规的计算方法彻底失效。为了不被这个"黑洞"吞噬，断裂力学家们发明了一种高明的战术：围道积分（Contour Integral）。我不直接去算黑洞中心的应力，而是在黑洞周围画一圈圈的"禁区"，通过计算跨过这些圈子的能量流，反推裂纹的破坏力。

断裂力学研究含裂纹结构的强度和行为，主要关注三个核心问题：

根据官方手册，要精准提取裂纹信息，Python 脚本需要完成以下三个关键动作。

`mdb.models[name].parts[name].engineeringFeatures.ContourIntegral(...)` 创建 Crack 对象。

你需要通过脚本指定裂纹的前缘（Crack Front）和延伸方向（Crack Tip/Bond Line）。

这就像是给台风定中心点。你得告诉 Abaqus 哪儿是刀锋，哪儿是刀背，算法才能围着这个点去转圈。

裂纹定义的关键要素：

1. 裂纹前缘（Crack Front）：

• 裂纹尖端的几何线

• 可以是边、节点或面

• 对于 3D 裂纹，是一条线

• 对于 2D 裂纹，是一个点

2. 裂纹尖端（Crack Tip）：

• 裂纹前缘的端点

• 用于定义局部坐标系

• 通常与裂纹前缘重合

3. 扩展方向（Extension Direction）：

• 裂纹预计扩展的方向

• 定义局部坐标系的 x 方向

• 影响 J 积分的计算

4. 裂纹面法向（Crack Normal）：

• 垂直于裂纹面的方向

• 定义局部坐标系的 z 方向

• 用于区分裂纹上下表面

3D 裂纹的定义：

def define_3d_surface_crack(model_name, part_name):    """    定义 3D 表面裂纹    场景：半椭圆表面裂纹    """    model = mdb.models[model_name]    part = model.parts[part_name]    part.ContourIntegral(        name='Surface-Crack',        crackFront=part.sets['Crack-Front-3D'],  # 裂纹前沿线        crackTip=part.sets['Crack-Tip-3D'],       # 裂纹最深点        extensionDirection=(0.0, 1.0, 0.0),       # 向内部扩展        symmetric=OFF,        crackNormal=(0.0, 0.0, 1.0)    )    print("3D 表面裂纹已定义")# 使用示例define_3d_surface_crack('3D-Fracture', 'Component')

通过参数 `numberOfContours` 设置围道层数。J 积分（能量释放率）在理论上与路径无关。

为了保险，我们通常会在裂尖外面画 5-10 层"圈"。如果第 3 层到第 10 层的计算结果都差不多，那就说明你的结果"收敛"了，是可信的。如果每层跳动很大，说明你的网格在"禁区"附近画得太糙了。

围道积分的收敛性分析：

理想情况下，J 积分与积分路径无关：

实际计算中，由于：

• 网格离散误差

• 数值积分误差

• 边界条件影响

不同围道的 J 积分会有差异。收敛性判据：

围道层数的选择：

围道位置的考虑：

• 第 1-2 层围道：太靠近裂尖，受奇异性影响大，通常不用

• 第 3-5 层围道：开始收敛，可用于初步判断

• 第 5-10 层围道：通常已收敛，推荐使用

• 外层围道：远离裂尖，可能受边界影响

在 Singularity 选项卡中设置中节点偏移（通常设为 0.25）。

为了在有限元网格里模拟出那种"应力无限大"的趋势，我们会人为地把裂尖单元的中节点往中心拨动一下。

这就像是给平淡的网格加了一点"重力加速度"，让原本线性的变化在裂尖处瞬间变陡。这 0.25 的偏移，就是断裂力学精度的灵魂所在。

偏移量的选择：

裂尖网格的最佳实践：

1. 辐射状网格（Spider Web）：

• 以裂尖为中心，向外辐射

• 第一层单元尺寸：a/10 ~ a/20（a为裂纹尺寸）

• 单元类型：二次单元（CPE8R, C3D20R）

2. 单元形状：

• 2D：三角形或四边形

• 3D：楔形或六面体

• 避免使用四面体（精度差）

3. 网格密度：

• 裂尖区域：极细

• 过渡区域：逐渐变粗

• 远场区域：正常密度

在 Python 中定义一个围道积分输出，你需要精准控制 engineeringFeatures：

from abaqus import *from abaqusConstants import *# 1. 在零件上定义裂纹特征# crackFront: 裂纹前缘的边或节点# crackTip: 裂纹尖端节点# extensionDirection: 裂纹预计的生长方向坐标p = mdb.models['Model-1'].parts['Plate']p.ContourIntegral(    name='Crack-1',     crackFront=p.sets['Crack-Edge'],     crackTip=p.sets['Crack-Tip-Node'],     extensionDirection=(0.0, 1.0, 0.0),     symmetric=OFF,     crackNormal=(0.0, 0.0, 1.0))# 2. 在分析步中请求输出（输出 J 积分和应力强度因子 K）mdb.models['Model-1'].fieldOutputRequests['F-Output-1'].setValues(    contourIntegral='Crack-1',     numberOfContours=5,     contourType=J_INTEGRAL)

别忘了在 Mesh 模块里，裂尖必须采用特殊的"蜘蛛网"辐射状网格（Ring of triangles/wedges）。

如果你用的是 XFEM（扩展有限元），那就更爽了，不需要画蜘蛛网网格，裂纹可以直接横跨单元生长！

完整的断裂分析流程：

from abaqus import mdbfrom abaqusConstants import *def setup_perfect_fracture_analysis(model_name):    model = mdb.models[model_name]    # --- 步骤 1：材料属性（力学 + 损伤） ---    mat = model.Material(name='Fracture-Material')    mat.Elastic(table=((210000.0, 0.3),))    mat.Plastic(table=((355.0, 0.0), (500.0, 0.2)))    # 定义延性损伤    mat.DamageInitiation(type=DUCTILE, table=((0.1, 2.0, 0.0),))     # --- 步骤 2：裂纹定义（移至 Interaction 模块） ---    # 假设已经在 Assembly 中创建了实例 'PART-1-1'    # 注意：crackFront 和 crackTip 应该是 Assembly 级别的集合    a = model.rootAssembly    model.ContourIntegral(        name='Main-Crack',        crackFront=a.instances['PART-1-1'].sets['Crack-Front'],        crackTip=a.instances['PART-1-1'].sets['Crack-Tip'],        extensionDirectionMethod=Q_VECTORS,        qVectors=(( (0.0, 0.0, 0.0), (1.0, 0.0, 0.0) ),),        symmetric=OFF    )    # --- 步骤 3：分析步 ---    model.StaticStep(name='Loading-Step', previous='Initial', nlgeom=ON)    # --- 步骤 4：修正输出请求（使用历史输出） ---    # 围道积分数据（J, K）必须通过历史输出请求获取    model.HistoryOutputRequest(        name='H-Output-J-Integral',        createStepName='Loading-Step',        contourIntegral='Main-Crack',        numberOfContours=5,        contourType=J_INTEGRAL # 可选 J_INTEGRAL, K_FACTORS, T_STRESS    )    print("断裂分析全流程已通过质检！")    return model

J 积分与应力强度因子的转换：

import mathdef calculate_stress_intensity_v2(j_integral, e_modulus, poisson_ratio, mode='I', state='strain'):    """    从 J 积分精确计算应力强度因子 K    参数:        j_integral: J 积分值 (N/mm)        e_modulus: 弹性模量 (MPa)        poisson_ratio: 泊松比        mode: 断裂模式 (目前支持 'I')        state: 'strain' (平面应变) 或 'stress' (平面应力)    返回:        K_m: 应力强度因子 (MPa·√m)    """    # 1. 根据约束状态确定有效弹性模量 E'    if state == 'strain':        # 平面应变        e_prime = e_modulus / (1 - poisson_ratio**2)    else:        # 平面应力        e_prime = e_modulus    # 2. 计算 K (单位: MPa·√mm)    # 公式: J = K² / E' -> K = sqrt(J * E')    k_mm = math.sqrt(j_integral * e_prime)    # 3. 单位换算: MPa·√mm -> MPa·√m (除以 sqrt(1000))    k_m = k_mm / math.sqrt(1000.0)    print(f"--- 质检计算结果 ---")    print(f"有效弹性模量 E': {e_prime:.2f} MPa")    print(f"应力强度因子 K: {k_m:.2f} MPa·√m (已完成单位转换)")    return k_m# 使用示例 (以典型钢材为例)# K_I = calculate_stress_intensity_v2(0.5, 210000, 0.3, state='strain')

在 Abaqus 脚本建模中，ContourIntegral 是我们与物理极限的一次握手。

最佳实践建议：

面对无法直视的"无穷大（奇异性）"，聪明的做法不是硬碰硬，而是绕道而行，通过"围观"来破解真相。掌握了围道积分的 API，你就掌握了量化"毁灭"的标尺。

围道积分的局限性：

• 仅适用于静态裂纹（不扩展）

• 需要预先知道裂纹位置

• 对网格质量要求极高

• 裂纹扩展需要重新网格化

👉互动话题：你在做断裂仿真时，遇到过 J 积分不收敛的情况吗？是网格太少，还是裂尖方向定反了？在评论区留下你的"排雷记录"，下一期我们聊聊如何让裂纹真正的"动起来"——XFEM 动态生长！

公众号	知识星球	扫码关注 科趣范，让科研更简单
小红书	抖音