近场动力学(PD)方法[1]由于采用空间积分方程来描述物质的运动平衡方程，而非传统的偏微分方程，因此特别适用于裂纹、断裂等非连续问题的求解．近场动力学方法获得迅猛发展，目前在材料本构模型、高精度数值求解方法及极区、岩土、热力、多物理场耦合等相关领域取得了一系列突破性的进展[2-13]．由于近场动力学方法采用非局部思想将结构进行均布空间离散，离散的物质点不只和相邻物质点发生作用，还会影响一定邻域范围内其他物质点的运动，导致其计算量要远大于传统的有限元方法(FEM)．陆续有研究者开展近场动力学和有限元耦合方法研究，致力于在发挥近场动力学求解非连续问题优越性的同时尽可能提升计算效率．文献[14]最早开始近场动力学和有限元耦合作用的研究，通过有限元中的杆单元类比近场动力学物质点对，实现了计算域整体刚度矩阵的表达，开展了钢球冲击韧性铝板和弹体破穿铁板等算例的分析，并且和Abaqus计算结果进行对比，结果显示了较好的符合度．文献[15]在此基础上将重合区域替换成界面单元来实现近场动力学和有限元的耦合，并且进一步对比了界面单元内等效节点载荷不同分配方案的结果差异．文献[16-17]提出了一种混合函数方法，实现了两种数值方法的耦合，结合近场动力学应变能表达函数和有限元应变能函数特征，提出了耦合区域的本构模型，较好地实现了局域和非局域数值计算模型的耦合．文献[18-20]提出了一种同一空间离散分布特征的近场动力学和有限元耦合模型，把近场动力学物质点和有限元节点统一起来，实现计算域的均匀离散，进而构建统一的刚度矩阵，实现不同计算域位移信息的同步求解．文献[21-22]采用同样的思路开展了二维裂纹扩张等准静态问题的求解，发现控制每一时间步下的键断裂数量，能在保证计算精度的同时显著提升计算效率．文献[23]也沿用了同样的思路，提出一种基于隐式格式的近场动力学和有限元耦合方法，并给出了非线性动态问题中的近场动力学方法和有限元法的等效增量方程，结合对Kalthoff-Winkler冲击问题研究，验证了数值预报方法的有效性．文献[24-27]同样采用全域空间均布离散、构建整体刚度矩阵求解的思路，开展了近场动力学和有限元方法耦合问题研究，显示基于该种思路的耦合方法具有广泛的应用性．文献[28]提出了一种通过力耦合的近场动力学与有限元混合建模的隐式分析方法，利用杆单元连接近场动力学子域与有限元子域，给出了杆单元的劲度表达式，并开展了含初始裂缝三点弯曲问题研究，验证了数值模型的有效性．综上可以发现文献[18]采用的全域空间均布离散的耦合思路得到了广泛应用，本研究也沿用这一思路，考虑到其方法及现有应用须要构建统一刚度矩阵方程，不仅实施步骤相对繁琐，且近场域的非局部特征决定了物质点在刚度阵表达中的非零元素数量较多，不具有带状稀疏性，容易给计算机内存带来较大负荷．本研究在此基础上提出了一种新的耦合交界面模型，针对不同计算域可以采用不同的求解器计算，通过设置交界面单元和虚拟边界层，连接近场域和有限元域实现混合建模，不用构建统一的刚度矩阵方程；并且进一步验证了交界区域可以通过形函数插值技术，实现由近场域细网格向有限元域粗网格的几何过渡；最后通过二维带预制裂纹板的拉伸和三维柔性杆的冲击等基础算例，验证了本研究提出的混合建模策略的有效性．1 耦合模型方法1.1　近场动力学方法近场动力学方法将计算域离散成有限个带有一定体积和质量的物质点．不同于传统的局域计算方法，物质点不仅和毗邻的粒子发生作用，还会通过键与一定邻域范围内(R=δ)的其他粒子形成物质点对，相互作用．为了更好地描述物质点之间的本构关系，将物质点和邻域物质点之间的相对位置定义为ξ=x'-x，相对位移矢量定义为η=u(x',t)-u(x,t)，u为该物质点在t时刻的位矢．基于虚功原理，推导可得近场动力学方法框架下物质点的运动平衡方程[7]ρu¨(x,t)=∫Hxf̲[x,t]x'-x-f̲[x',t]x-x'dVx'+b(x,t),(1)式中：ρ为材料密度；Hx为近场邻域空间范围；Vx'为近场邻域体积；b为外力密度；f为物质点对之间的本构力，可通过物质点对之间的键常数与键伸长率来表达求解．有f(η,ξ)=(ξ+η)/ ξ+ηcs(ξ)μ(x,ξ)，(2)式中：c为键常数；s(ξ)为键断裂常数，s(ξ)=ξ+η-ξ/ ξ;(3)μ(x,ξ)为键断裂破坏的判断函数，当键的伸长率超过其极限阈值s0时键就会断裂，意味着粒子点对失去力的作用，μ(x,ξ)=1    (s(ξ)s0),0    (其他).(4)点的局部损伤程度表示与该物质点相连的所有键断裂的程度，用局部损伤函数φ表示，φ(x,t)=1-∫Hxμ(x,t,ξ)dVx'/∫HxdVx'．(5)1.2　有限元动力方程有限元动力方程可以写为[29-30]Ma¨+Ca˙+Fint=Fext，(6)式中：M为结构刚度矩阵；C为阻尼矩阵；a¨和a˙为有限元节点的加速度和速度；Fint和Fext为结构所受到的内力和外力．又有Fint=∫VBTσdV;(7)Fext=∫VNTρbdV+∫ANTt¯dA,(8)式中：B为应变位移矩阵；σ为内应力；N为节点形函数．当进行时域积分时，往往采用隐式和显式两种求解格式，考虑到本研究在近场动力学求解格式中选择了显式格式，这里同样采用显式时域积分的中心差分求解格式求解．已知0，t1，t2，…，tn时刻的位移、速度和加速度，则n+1时刻的位移an+1和n+1/2时刻的速度a˙n+1/2可表示为：an+1=an+Δtn+1/2a˙n+1/2；(9)a˙n+1/2=a˙n-1/2+Δtna¨n．(10)经过转换可以得到：a˙n+1/2=M/Δt-C/2M/Δt+C/2a˙n-1/2+Fext-FintΔtM/Δt+C/2；(11)an+1=an+1+a˙n+1/2Δt．(12)1.3　近场动力学和有限元高效耦合求解策略提出一种高效的近场动力学和有限元耦合策略如图1所示，其中物质点尺寸为Δx，有限元网格尺寸为ΔL．计算域被划分成四块区域，纯近场动力学域如图1中绿色点所示；近场域和有限元域的交界面区域如图1中黄色点所示，其有两种属性，既是近场动力学物质点，又是有限元节点；近场域的虚拟边界层区域如图1中粉红色点物质点所示，其有两种属性，既是近场动力学虚拟边界层上的物质点，又是有限元节点；纯有限元区域如图1中黑色节点所示．其中近场动力学和有限元耦合方法主要通过交界面和虚拟边界层上的节点信息交互传递实现，其具体传递思路如下．10.13245/j.hust.241221.F001图1近场动力学和有限元耦合模型a．近场动力学域到有限元域的运动信息传递近场域和有限元域的交界面区域上的黄色点有两种属性，既是近场动力学物质点，又是有限元节点，两种点共享同一空间位置．通过将交界面上近场动力学物质点的位置和速度信息传递到同一空间位置的有限元节点上，实现了近场域到有限元域运动信息的传递．当其被划分成近场域或有限元域时，该点的运动须要遵循对应计算域应满足的运动平衡方程．b．有限元域到近场动力学域的运动信息传递近场域的虚拟边界层上的粉红色点有两种属性，既是近场动力学虚拟边界层上的物质点，又是有限元节点，两种点共享同一空间位置，通过将该种属性上的有限元节点的位置信息传递到同一空间位置上的近场动力学虚拟边界层的物质点上，实现了有限元域到近场域的运动信息的传递．耦合方法的主要流程步骤如下．步骤1 结构几何、材料、初始运动参数输入，不同计算域网格模型建立；步骤2 近场域内邻域粒子搜索，表面因子和体积修正；步骤3 时域循环，由有限元耦合区域网格节点位移场确定虚拟边界粒子位置信息；步骤4 近场域内物质点键力计算，物质点位置信息更新；步骤5 物质点速度场、位移场传递到交界面有限元节点，有限元域内节点结构动力计算；步骤6 更新有限元域内节点运动信息，判断是否达到最大计算步长，是则结束，否则重复步骤3～6．2 算例分析2.1　柔性杆的冲击算例一为两根相同的柔性杆冲击问题模拟，如图2所示．杆的长度L=0.05 m，宽度W=0.01 m，高度H=0.01 m；左右杆的冲击速度为±10 m/s；材料的弹性模量E=75 GPa，泊松比υ=0.25，密度ρ=10.13245/j.hust.241221.F002图2柔性杆的冲击示意图2 700 kg/m3．采用纯近场动力学方法和近场动力学-有限元方法(PD-FEM)耦合两种方法进行模拟．左右端1/5的杆长采用有限元域方法模拟，中间区域采用近场动力学粒子方法模拟．不同计算域采用统一的网格尺度，为了更好地对比验证，本研究近场域的物质点半径保持和文献[7]一致，即Δl=Δx=0.1 mm．时间步长选择为Δt=9.318×10-8 s．图3为分别采用近场动力学、近场动力学-有限元耦合方法计算的杆结构X方向位移场对比．图4为左右杆节点的轴向位移场的时域变化情况，并与文献[7]中有限元分析(FEA)计算结果进行了对比．从图3可以很直观地发现纯近场动力学方法计算的位移场和耦合模型计算结果十分符合，特别是在耦合边界区域位移场由近场动力学计算域到有限元计算域过渡十分顺滑，表明不同计算域位移协调性良好．图4的节点位移曲线也显示了用本研究提出的耦合方法得到的左右杆的中心节点位移变化和纯近场动力学方法及有限元分析方法的计算结果十分相似．首先在初始速度作用下左右杆会相向运动，而后在第56～141计算步时域范围内左右杆相互接触，给彼此一个与运动方向相反的作用力，中心节点沿X方向位移几乎保持不变，最后超过141计算步后左右杆件会朝着相反的方向分离运动．10.13245/j.hust.241221.F003图3近场动力学方法与近场动力学-有限元耦合方法结果对比(t=9.318×10-7s)10.13245/j.hust.241221.F004图4不同计算方法得到的节点(±0.012 5，0.0，0.0)轴向位移时历变化对比2.2　带预制裂纹板的拉伸算例二为含I型裂纹的二维平板在拉力作用下的破坏．板的尺寸为L1×L2=50 mm×50 mm．板的中间有一条沿着x方向分布的长度为a1=10 mm的预制裂纹．板的上下端受到拉力外载作用，以V=±20 m/s的速度匀速运动．材料的弹性模量为E=192 GPa，泊松比为υ=0.33，密度ρ=8 000 kg/m3．同样采用纯近场动力学方法和近场动力学-有限元耦合两种方法进行模拟，其中上下端20%板宽采用有限元域方法模拟，中间区域采用近场动力学粒子方法模拟．不同计算域采用统一的网格尺度，即Δl=Δx=0.001 m，时间步长选择为Δt=1.336 7×10-8s．图5和图6分别给出了采用近场动力学方法和近场动力学-有限元耦合方法计算的拉力外载下板的位移云图和损伤云图，可以看到不同时刻下耦合方法和近场动力学方法计算的位移云图十分符合，另外耦合方法计算位移在不同模型区域的过渡十分光滑．从板的损伤云图可以发现近场动力学方法和耦合方法计算的裂纹扩展过程一致，在拉伸外力作用下，裂纹沿着尖端在X方向不断扩张，这表明采用耦合方法能保留近场动力学方法在模拟裂纹扩展等非连续特征上的优势．10.13245/j.hust.241221.F005图5近场动力学方法和近场动力学-有限元耦合计算方法Y方向位移云图对比(t=1.17×10-5 s)10.13245/j.hust.241221.F006图6近场动力学方法和近场动力学-有限元耦合方法计算的板的损伤云图(t=1.17×10-5 s)为了更直观地反映计算域局部区域的结果差异，跟踪记录了上下两个粒子点(0，±0.01)全时域的位移情况，如图7所示．可以看到：由于加载区域和记录的坐标点存在一定距离，应力波传递到坐标点区域需要一定时间，须经大约224次计算步后上下粒子点才沿着各自的加载方向发生位移，因此图7在前半段的粒子点位移一直在0刻度线附近．记录的上坐标点和下坐标点朝着反向运动，这是因为板块上下边界施加反向的速度导致．从整体上看，两种模型计算结果差异较小，其中最大误差为2.6%，平均误差为2.2%，两种模型计算结果显示出整体时域范围都比较符合的特征．这也进一步支撑了本研究提出的近场动力学-有限元耦合求解策略预报的有效性．10.13245/j.hust.241221.F007图7近场动力学方法和PD-FEM耦合方法计算点的位移本研究提出的耦合思路须在交界面实现由近场动力学物质点运动信息到有限元节点运动信息的传递，Galvanetto等[14]采用统一的空间离散域特征，即有限元域和近场域采用相同的网格尺度，这可能导致有限元网格划分过于精细，增加了计算时间消耗．图8为有限元计算域采用不同网格尺寸预报的板Y方向的位移，可以看到有限元网格尺度变化对计算结果影响不大，这是因为近场域的网格精度依赖性要大于有限元域，在实际计算求解特别是大型工程问题的求解中可以适当调大有限元网格单元大小以减少计算量．计算所用的电脑处理器型号为Intel(R) Core(TM) i7-10875H CPU@2.30 GHz．通过对比发现，纯近场动力学方法的计算耗时为866 s，而同样空间分布离散特征的近场动力学-有限元耦合算法的计算耗时为654 s，计算效率提升24%，显示了耦合算法在计算速度上的优势；另外，当分别采用2倍和4倍粒子半径的有限元网格尺寸单元时，计算耗时分别为491 s和383 s，相比于1倍粒子半径的计算耗时得到大幅度缩减，且计算结果彼此符合度较好，表明本研究采用的求解策略实现由近场域细网格向有限元域粗网格的几何过渡，能减少耦合算法中冗余网格的计算耗时，进一步提升计算效率，为近场动力学-有限元耦合提供了新的思路，显示出较好的工程应用潜力．10.13245/j.hust.241221.F008图8近场动力学-有限元耦合方法预报板Y方向位移云图3 结论针对近场动力学方法计算存在的效率问题，提出了一种新的高效的近场动力学和有限元耦合策略，主要得到以下结论．a．本研究利用建立的近场动力学-有限元耦合模型开展了柔性杆冲击算例分析，计算结果和纯近场动力学方法预报结果和有限元分析结果符合度较好，且耦合边界区域位移场由近场动力学计算域到有限元计算域过渡十分顺滑，显示了良好的位移协调性，验证了本研究提出的耦合求解策略准确可靠．b．当采用耦合模型预报二维预制裂纹板裂纹扩展模拟时，裂纹沿着尖端在X方向不断扩张，与纯近场动力学方法预报结果一致，表明采用耦合方法能保留近场动力学方法在模拟裂纹扩展等非连续特征上的优势．c．当有限元网格尺寸和近场动力学粒子直径相等时，耦合模型计算效率提升24%，证明了耦合算法在计算速度上的优势；当分别采用2倍和4倍粒子半径的有限元网格尺寸单元时，计算效率得到进一步提升，且计算结果彼此符合度较好，表明本研究采用的由近场域细网格向有限元域粗网格几何过渡的求解策略，能减少耦合算法中冗余网格的计算耗时，进一步提升计算效率．本研究中近场动力学-有限元耦合模型不能实现近场动力学域粒子和有限元节点的实时转换，因此近场动力学域的范围须要依据计算问题中的裂纹区域和裂纹扩展路径进行选取．