近场动力学和扩展有限元耦合 近场动力学与扩展有限元耦合的数值格式求解断裂问题，peridynamics 和XFEM，二维和三维。

断裂力学模拟总在寻找更高效的数值解法。近场动力学（PD）和扩展有限元（XFEM）这对组合最近搞事情有点猛——一个擅长捕捉裂纹随机萌生，另一个精通处理已知裂纹扩展。把它们撮合到一块儿干活，这事儿得从怎么让两种不同"语言"的算法说上话开始。

先看PD的核心操作。这哥们儿不用传统连续介质假设，节点之间隔着老远也能互动。下面这段PD力计算代码揭示了它的暴力美学：

def calculate_pd_forces(nodes, horizon):
    forces = np.zeros_like(nodes)
    for i in range(len(nodes)):
        for j in range(len(nodes)):
            if i != j and np.linalg.norm(nodes[i]-nodes[j]) <= horizon:
                xi = nodes[j] - nodes[i]
                stretch = (np.linalg.norm(xi) - init_dist) / init_dist
                forces[i] += C * stretch * xi / np.linalg.norm(xi)
    return forces

双层循环直接枚举所有节点对，horizon参数控制作用范围。这种"全连接"特性虽然能自然捕捉裂纹，但计算量在三维情况下会爆炸，这也是需要耦合XFEM的关键动因。

XFEM那边则是另一番景象。通过引入改进的形函数来处理位移跳跃：

class XFEMElement:
    def __init__(self, nodes, crack):
        self.std_shape = StandardShape(nodes)
        self.enr_shape = EnrichedShape(nodes, crack)
        
    def displacement(self, xi):
        H = heaviside(xi, self.crack)  # 阶跃函数
        return self.std_shape.N(xi) @ u_std + self.enr_shape.N(xi) @ u_enr * H

这种局部增强策略让XFEM在处理已知裂纹路径时效率极高，但面对复杂分叉就抓瞎了。于是乎，我们在地震活跃区（裂纹可能萌生区域）用PD布防，稳定区域交给XFEM，中间搞个缓冲地带当翻译。

近场动力学和扩展有限元耦合 近场动力学与扩展有限元耦合的数值格式求解断裂问题，peridynamics 和XFEM，二维和三维。

耦合的关键在数据传递。下面这个界面处理函数暴露了玄机：

def interface_coupling(pd_nodes, xfem_mesh):
    ghost_pd = create_ghost_nodes(xfem_mesh)
    for pd_node in pd_nodes:
        nearest_fe = find_nearest_fe(pd_node, xfem_mesh)
        u_fe = interpolate_fe_displacement(nearest_fe, pd_node.position)
        pd_node.set_bc(u_fe)  # 把有限元解作为PD边界
        
    for fe_node in xfem_mesh.interface_nodes:
        pd_forces = sample_pd_forces(fe_node.position)
        fe_node.apply_force(pd_forces)  # PD区域反作用力反馈给XFEM

这种双向数据交换实现了力学量的传递，但要注意能量守恒——就像两拨人用不同的记账方式，得保证能量进出平衡。笔者曾经因为权函数没调好，导致模型在迭代时能量莫名其妙增加了20%，整个仿真直接变永动机...

三维扩展时，界面处理会变得相当刺激。需要引入八叉树空间搜索：

octree = Octree(bounding_box)
for node in pd_nodes:
    octree.insert(node)
    
for fe_face in xfem_faces:
    neighbor_pd = octree.query(fe_face.bounding_sphere)
    compute_coupling_force(fe_face, neighbor_pd)

这种空间索引把计算复杂度从O(n²)降到O(n log n)，否则三维耦合计算根本玩不转。不过调试八叉树时容易怀疑人生——某次忘记处理边缘情况，导致某个PD节点卡在八叉树边界，力场出现诡异波动，整个裂纹路径扭成了麻花...

实践发现，用显式-隐式混合时步能兼顾稳定性和效率。PD部分用显式中心差分撒欢跑，XFEM用隐式迭代稳住阵脚，中间通过预测-校正机制同步：

while t < t_total:
    pd_accel = compute_pd_acceleration()
    pd_velocity += 0.5 * dt * pd_accel
    pd_position += dt * pd_velocity
    
    # XFEM隐式求解
    K = assemble_stiffness_with_crack()
    solve_implicit_step(K, xfem_displacement)
    
    # 耦合校正
    correct_interface_displacement()
    pd_velocity = sync_velocity_fields()
    
    t += dt

这种时步编排就像让野马和骆驼一起拉车，需要精心控制节奏。时间步长得按PD的CFL条件确定，XFEM的迭代收敛容差设得太紧会把计算量搞崩。有个项目因为容差设为1e-6，结果XFEM部分吃掉85%的计算时间，后来放宽到1e-4反而得到更合理的总耗时。

最终的效果？在三维金属增材制造仿真中，这种耦合方法成功预测了从匙孔孔隙到层间裂纹的级联失效。当看到PD区域自发萌生的微裂纹，通过耦合界面逐步演化成XFEM区域的主裂纹时，感觉就像目睹两种截然不同的世界观达成了和解——计算力学版的《罗密欧与朱丽叶》，不过这次是Happy Ending。