本发明属于计算力学领域，提供了一种结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，该方法将时间间断的思想引入态基近场动力学理论，有效提高了近场动力学显式动力分析的精度和准确预测结构断裂破坏的能力。本发明采用时间间断显式时程积分格式可以有效控制传统时程积分方法带来的虚假数值振荡，采用非常规态基近场动力学模型简便全面地描述了材料在冲击荷载下的复杂力学行为，并通过多种损伤断裂准则有效表征了结构的冲击断裂破坏模式。此外，本发明还采用快速邻域搜索算法构建物质点邻域并更新接触邻域，提升了计算效率。本发明所提出的方法作为一种新的数值求解格式，可以通过简单修改原计算程序实现，降低了数值实施复杂度。
1.一种结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，其特征在于，步骤如下：首先，给出非常规态基近场动力学求解的基本格式；基于非常规态基近场动力学框架，控制方程以积分‑微分方程的形式给出，其中，t表示时间，ρ是密度，x是参考构型中的物质点，u是位移矢量，是加速度，T是力状态，b是给定的外力密度，Vx′是物质点x’的体积，Hx是半径为δ的物质点x的邻域；定义ξ＝x’‑x是参考构型中的相对位置，y和y’是当前构型中物质点的位置，η＝u’‑u是物质点x和x’之间的相对位移；在非常规态基近场动力学中力状态T与变形状态Y＝y’‑y＝η+ξ并不平行；
非常规态基近场动力学将经典材料本构模型引入近场动力学框架，它基于经典变形梯度F的非局部表征，其中，K是形状张量， 表示两矢量并乘，
这里物质点之间相互作用的程度由标量函数ω描述；在连续介质力学中，速度梯度从变形梯度导出，得到非局部变形梯度F和非局部速度梯度L后，结合经典材料本构模型，得到一系列应力应变张量；这样，非常规态基近场动力学力状态为，其中，Px表示物质点x的PK1应力，Kx表示物质点x的形状张量；
该方程使用数值方法求解；系统在空间上离散为Nx个物质点，xp的邻域有Np个物质点；
时域I＝[0，T]被均匀划分为Nt个区间，离散时刻为tn＝nT/Nt，其中n＝0，1，…，Nt；物质点xp在tn时刻的控制方程(6)写为，这里下标p和q表示点的编号；至此，非常规态基近场动力学的求解框架就建立起来；
接着，将经典材料本构关系引入非常规态基近场动力学求解框架；由式(2)求得变形梯度F后，通过计算出格林应变张量εij；其中δij为克罗内克尔符号，这里的下标是张量的指标；在得到时间步n的格林应变张量εij后，已知物理量是上一步中的格林应变张量 塑性应变张量内变量 本步应变增量 这里上一步中的物理量用下标n‑1标出，本步骤中的物理量在下标n中省略；塑性问题的数值计算用返回映射算法，其由两个步骤组成，弹性试验步假设材料具有弹性，试应力 Cijkl表示弹性张量，表示上一步应力张量，Δεkl表示本步应变增量；此时的应力偏离该步骤的真实屈服面，然后使用径向返回算法进行塑性校正，将试应力拉回到屈服面上；对于J2流动理论的特殊情况，只有一个内部变量 屈服条件是其中，sij是偏应力， 是等效塑性应变，σy是当前屈服应力；径向返回算法一般使用循环迭代判断的方法实现，第k个迭代步的应力和塑性应变的增量为其中，nij是屈服面的法向方向，且
其中，G是剪切模量，E是塑性模量，f是第k步的屈服函数值，当|f|＜TOL时迭代结束，其中TOL为给定的迭代收敛容差；
根据Johnson和Cook提出的有关应变硬化、应变率和热软化的流动应力模型，则应变率效应必须包含在冲击和爆炸问题的本构模型中，其中，A、B、C、m1和m2是Johnson‑Cook模型的材料常数； 是塑性应变率， 是有效塑性应变率，Tr和Tm是室温和熔化温度，T是当前温度；Johnson‑Cook模型还定义了物质点随时间的累积损伤，其中， 为等效塑性应变增量，当D＝1时物质点发生断裂； 是断裂的等效应变，其中，σm是三个法向应力的平均值，是Mises等效应力，D1、D2、D3、D4和D5是Johnson‑Cook损伤的材料常数；通过损伤来折减应力，当损伤参数D达到1时，发生断裂；断裂发生后，偏应力置零；如果球应力为正，则再将球应力置零；
最后，导出时间间断态基近场动力学的基本格式；在时间间断伽辽金方法中，允许未知场函数在相邻时间间隔之间间断；函数w(t)在时间域的阶跃定义为，其中， 是间断算子符号；对
于任意的时间间隔In，其初时刻和末时刻分别为tn和tn+1，位移u和速度v分别基于三次Hermite函数和线性函数进行插值，其中，vn和vn+1是本时间步tn和tn+1时刻的速度，un和un+1是本时间步tn和tn+1时刻的位移，且时间步长Δt＝tn+1‑tn，I是单位矩阵；类似地，假设外力密度b在每个时间间隔内线性变化，内力密度f通过Hermite函数独立插值，其中，bn和bn+1是本时间步tn和tn+1时刻的外力密度，fn和fn+1是本时间步tn和tn+1时刻的内力密度， 和 分别表示内力密度在本时间步tn和tn+1处的时间导数；
虽然物质点的位移和速度是独立插值的，但它们必须满足如下约束条件，运动方程，约束条件，连同时间间隔In内位移和速度的间断方程一起构成积分弱形式，从δvn、δvn+1、δfn和δfn+1的独立变分得到SBPD‑TD的基本公式，即，其中，ρ是密度， 是上一步计算出的末时刻位移， 是上一步计算出的末时刻速度，是上一步计算出的末时刻内力，且以上的求解格式是一种新的逐步积分法，用以替代求解动量方程(7)常用的中心差分法，从而控制由中心差分法带来的应力解中的数值振荡；
根据上述理论推导得出的时间间断态基近场动力学的基本格式，采用显式求解策略，再结合经典材料本构计算方法即实现结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，具体步骤如下：步骤1：建立离散物质点模型，定义材料参数，根据离散的物质点位置坐标和预设的邻域大小判断每个物质点的邻域范围，建立物质点之间的相互作用关系；
步骤2：定义边界条件和时间步长，根据位移场和速度场利用态基近场动力学方法计算结构的内力和内力的时间导数；引入经典的塑性本构及损伤演化法则进行非线性分析，更新应力和塑性应变；
步骤3：通过实施交错迭代求解策略，利用式(26)‑(29)进行显式时程积分，更新位移和速度，当迭代收敛时获得本步的位移解和速度解，返回步骤2进行下一时间步计算。
2.根据权利要求1所述的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，其特征在于，在步骤1中，为了更加精确地贴合复杂结构的表面，使用有限单元法中的六面体网格进行离散，单元的形心作为物质点坐标，单元的体积作为物质点体积；作为一种非局部模型，近场动力学的物质点与其邻域内所有的物质点都具有相互作用；为了减少计算规模，基于胞元构建邻域；物质点x的视界定义为以δ为特征的范围内所有物质点；划分胞元后，只需搜索x所在胞元及其周围8个胞元中的物质点即得到物质点的邻域；对于三维情况，需要搜索
3.根据权利要求1所述的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，其特征在于，在步骤2中，为了抑制非常规态基近场动力学中固有的零能模式，引入一种基于键的变形梯度修正方法，定义变形状态的修正为zpq＝Ypq‑Fpξpq (31)
其中g是正常数，c0是键基近场动力学中的微模量，ω0定义为
4.根据权利要求1所述的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，其特征在于，在步骤2中，使用径向返回算法更新应力和塑性应变；J2流动理论的屈服面在π平面内是一个圆，其法向沿圆的径向；采用关联塑性流动法则，塑性流动方向为，其中， 这里 表示试应力，
5.根据权利要求1所述的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，其特征在于，在步骤2中，对于非线性分析，本构模型用率形式表示；柯西应力张量σij的物质时间导数与Jaumann应力率 有关，其中，Ωij是旋率张量，它是速度梯度L的反对称部分，对于各向同性弹性模型，Jaumann应力率从变形张量率 导出，其中，Cijkl是弹性张量。
6.根据权利要求1所述的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，其特征在于，在步骤2中，引入两种损伤的形式；一种是基于键的损伤，一种是基于点的损伤；对于基于键的损伤，通过如下的三种方式判断键的断裂，(1)定义物质点p和q之间键的伸长率为
其中ξpq是物质点p和q在参考构型中的相对位置，ηpq是物质点p和q之间的相对位移；键的断裂状态μ由下式给出，其中，s0表示键的临界伸长率；
其中，Dc≤1是渐进性延性损伤材料参数，Dc＝1对应于键的突然断裂；Dξ是键的损伤参数，通过两端物质点的Johnson‑Cook损伤加权平均获得；
基于上述三种方法，在基于键的损伤模型中，近场动力学的标量影响函数ωpq需要修改为ωpq＝μωpq；此外，还需引入了一个标量函数 来描述物质点xp的失效程度，对于基于点的损伤，通过Johnson‑Cook模型定义物质点的损伤后，可以通过损伤来降低应力；当损伤参数D达到1时，发生断裂；断裂发生后，偏应力s设置为0；如果计算的球应力为正，则将其重置为零。
7.根据权利要求1所述的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，其特征在于，在步骤3中，式(26)‑(29)的显式时程积分方案使用不动点迭代的方法完成；具体实施过程使用伪代码的形式给出：(1)初始化变量，令循环变量k＝0；
(2.1)根据近场动力学方法，使用速度 和位移un、 分别计算内力密度的时间导数 和
(2.2)根据近场动力学方法，使用位移 计算内力密度(2.3)根据式(26)‑(29)，计算
(2.4)若 转进入步骤(3)，否则k＝k+1，返回步骤(2)继续循环；
8.根据权利要求1所述的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，其特征在于，结构冲击弹塑性断裂分析的时间间断态基近场动力学方法的具体实现过程将通过如下的伪代码形式展示：(1)建立离散的物质点模型，构建邻域，定义材料参数：密度ρ、弹性模量E、泊松比v、塑p性模量E 、Johnson‑Cook参数、初始屈服应力σy0、临界伸长率s0、临界等效塑性应变 渐进性延性损伤材料参数Dc，施加载荷并初始化速度位移等各变量；
(2.2.1)使用速度 和位移un、 分别计算内力密度的时间导数 和 以及内力
(2.2.3)若 进入步骤(2.3)；否则k＝k+1，返回步骤(2.2)继续循环；
(2.3)更新速度 位移 若n＜Nt，则n＝n+1，返回步骤(2)；否则时间步循环结束，进入步骤(3)；
结构冲击弹塑性断裂分析的时间间断态基近场动力学方法
技术领域
[0001] 本发明属于计算力学领域，具体涉及一种结构冲击弹塑性断裂分析的时间间断态基近场动力学方法。
背景技术
[0002] 冲击破坏问题广泛存在于工程实际中，如地震导致的山体滑坡、军事武器研制中子弹与靶体的相互作用、生活中电子产品的跌落、以及极地科考中破冰船与冰的相互作用。可以说工程事故中大概率存在冲击载荷和冲击动力学问题；而对于军用民用产品的研发，材料动态断裂力学行为的研究更是不可避免的。所以结构的冲击断裂力学行为分析在防灾减灾、武器研发、产品研制等方面发挥着重要的作用，对解决我国关键技术“卡脖子”难题和保障人民群众生命财产安全具有重要的意义。
[0003] 但是由于冲击破坏问题涉及复杂的非线性力学行为，例如脆/塑性断裂破坏、材料大变形、材料间接触等。对于冲击断裂破坏机理的研究，若通过实验手段不仅实施困难、代价高昂，而且很难同时考虑多种外部因素对断裂破坏行为的影响。数值仿真方法作为新兴高性能数值计算技术为深入探究结构的复杂断裂过程提供了一种可行的方案，其在选用合理的本构模型和有效的损伤法则后即可快速、准确地预测结构的动态断裂失效行为。因此，发展高精度、简便的材料弹塑性冲击断裂分析的数值模拟方法尤为重要。
[0004] 研究表明，结构和材料冲击断裂力学行为与其静态力学行为具有很大的不同。首先，材料和结构受到动载荷作用时，所产生的应力和变形将以波的形式传播开来，其变形和应力的不均匀性是不可忽略的，因而针对应力波的反射、透射、弥散现象的精确分析具有重要的意义。其次，在冲击爆炸等强动载荷作用下，结构中的材料将发生高速变形。根据材料的微观变形机制分析，材料对高速变形的抵御能力高于对缓慢变形的抵御能力，也就是材料在高速变形时通常具有明显的率相关性。再次，根据达朗贝尔原理，当动载荷超过极限静载荷时，结构的惯性力参加承受外载荷，抵抗变形，因而结构可以在短时间内承受比静极限载荷高得多的外载荷。最后，材料的断裂也会与结构的动力学行为相互影响，断裂动力学的问题在数学和物理上比断裂静力学要复杂得多，实际上是几种非线性现象相互耦合的问题。上述特点的存在，使得材料的冲击断裂行为分析目前仍然存在很大的挑战，具有重要的研究价值。
[0005] 目前材料冲击断裂力学行为的数值分析存在损伤断裂表征困难、计算精度不高的问题。在损伤断裂的表征方面，由于目前工程上常用的数值方法一般基于连续介质力学，在处理断裂等具有奇异性特征的问题时其空间微分无法正确定义。因此传统断裂力学一般将断裂行为归结为J积分、应力强度因子等参数来间接描述结构中的裂纹，且在动态裂纹扩展和弹塑性断裂行为的表征上仍然存在许多困难。而在计算精度方面，由于传统的时间离散格式一般基于差分法，所求得的应力解具有严重的数值振荡。这些虚假的数值振荡又进一步造成虚假的损伤，从而严重干扰正常损伤演化和裂纹扩展的计算。
[0006] 最近的研究表明，近场动力学模型是一种处理非局部损伤和裂纹扩展问题的有效方法，其通过建立空间上的积分方程表征材料的力学行为，从而解决了连续介质力学的微分方程在裂纹处不可求导的难题。近场动力学最早于2000年由Silling提出，并随后应用到了材料的裂纹扩展问题分析中。经过二十多年的快速发展，近场动力学因具有能准确裂纹奇异性、模型简洁、编程方便等优点而得到了越来越多学者的青睐，已经被应用于解决包括准静态断裂、动态断裂、塑性断裂、多相及多物理场断裂等在内的诸多具有奇异性和非局部性特征的工程问题。值得注意的是，近场动力学目前可以归结为两种模型，即键基近场动力学模型和态基近场动力学模型。前者由于具有泊松比的限制一般用于分析脆性材料的断裂力学行为；而后者由于消除了泊松比的限制且可以方便地引入传统弹塑性本构，可以用于材料弹塑性断裂力学行为的分析。因此，本发明将在空间上考虑使用态基近场动力学理论来研究弹塑性材料的冲击断裂力学行为。
[0007] 为了精确地模拟在受到冲击载荷时材料内应力波的传播过程，以准确预测裂纹的萌生与扩展，本发明在时间域上拟使用时间间断的思想。由于传统的时程积分方案一般采用差分方法，其速度和位移在时间域内是连续的。当受到冲击碰撞等存在不连续性特征的荷载时，传统时程积分方案不能表征响应的高梯度特征并会出现虚假的数值振荡。为了解决这个问题，时间间断的思想就被提出。它的主要特征是位移和速度在时间域内独立地插值，并允许它们在离散的时刻出现间断。这样时间间断方法就可以准确再现结构中的间断特征并抑制由传统积分方案所带来的虚假数值振荡。
[0008] 目前，时间间断方法的研究主要集中在结构中应力波的传播过程研究，关于冲击过程中应力波引起的结构断裂力学行为研究尚未充分开展，结构断裂行为和应力波传播耦合现象的数值模拟研究还不充分。因此，本发明将提出一种时间间断的态基近场动力学方法，控制传统数值方法所固有的数值振荡，以模拟结构在受到冲击载荷时发生的弹塑性断裂力学行为，为工程断裂问题的精确分析提供一种行之有效的途径。
发明内容
[0009] 本发明要解决的技术问题：本发明基于时间间断思想和态基近场动力学方法，并结合经典本构模型，针对结构冲击载荷下应力波引起的损伤断裂展开研究，创新性地提出了一种时间间断态基近场动力学方法(Time‑discontinuous state‑based peridynamics，简称SBPD‑TD)，其目的在于解决现有技术存在的以下问题：采用时间间断方法以克服传统时程积分格式难以准确捕捉应力解的高梯度特征和存在虚假数值振荡的缺点；采用近场动力学模型以克服传统连续介质力学框架下难以描述断裂行为中固有奇异性的不足，并降低数值实施的难度；克服基于传统有限单元法(FEM)模拟材料有限变形弹塑性断裂破坏时因发生网格畸变而造成数值精度严重降低的缺点；采用态基近场动力学方法弥补键基近场动力学方法所带来的泊松比固定且不易引入经典本构关系的不足；采用Johnson‑Cook模型计算材料的冲击响应并引入多种简便损伤表征，减少了传统断裂力学损伤表征的复杂度。
[0010] 本发明的技术方案：
[0011] 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法(SBPD‑TD)，
[0012] 为了准确捕捉冲击过程中结构应力的高梯度变化并抑制数值振荡，提高断裂分析的精度，本发明将基于态基近场动力学模型耦合时间间断方法分析冲击载荷下结构中应力波的传播特性并预测裂纹的萌生与扩展，提供了一种结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，具体步骤如下：
[0013] 首先，给出非常规态基近场动力学求解的基本格式。基于非常规态基近场动力学框架，控制方程以积分‑微分方程的形式给出，
[0014]
[0015] 其中，t表示时间，ρ是密度，x是参考构型中的物质点，u是位移矢量，ü是加速度，T是力状态，b是给定的外力密度，Vx′是物质点x′的体积，Hx是半径为δ的物质点x的邻域。图2′给出了态基近场动力学(SBPD)中的参考构型和当前构型，其中ξ＝x‑x是参考构型中的相′ ′
对位置，y和y是当前构型中物质点的位置，η＝u‑u是物质点x和x′之间的相对位移。值得注′
意的是，在非常规态基近场动力学中力状态T与变形状态Y＝y‑y＝η+ξ并不平行。
[0016] 非常规态基近场动力学的一个显著优势是它可以便捷地将经典材料本构模型引入近场动力学框架，它基于经典变形梯度F的非局部表征，
[0017]
[0018] 其中K是形状张量， 表示两矢量并乘，
[0019]
[0020] 这里物质点之间相互作用的程度由标量函数ω描述。在连续介质力学中，速度梯度可以从变形梯度导出，
[0021]
[0022] 得到非局部变形梯度F和非局部速度梯度L后，结合经典材料本构模型，可以得到一系列应力应变张量(如格林应变张量和PK1应力)。这样，非常规态基近场动力学力状态为，
[0023]
[0024] Px表示物质点x的PK1应力，Kx表示物质点x的形状张量。
[0025] 将式(5)代入式(1)，控制方程可以写作，
[0026]
[0027] 该方程可以使用无网格方法等数值方法求解。系统在空间上离散为Nx个物质点，xp的邻域有Np个物质点。时域I＝[0,T]通常被均匀划分为Nt个区间，离散时刻为tn＝nT/Nt，其中n＝0,1,…,Nt。物质点xp在tn时刻的控制方程(6)可写为
[0028]
[0029] 这里下标p和q表示点的编号。至此，非常规态基近场动力学的求解框架就建立起来。
[0030] 接着，将经典材料本构关系引入非常规态基近场动力学求解框架。由式(2)求得变形梯度F后，可通过
[0031]
[0032] 计算出格林应变张量εij，其中δij为克罗内克尔符号，这里的下标是张量的指标。在得到时间步n的格林应变张量εij后，已知物理量是上一步中的格林应变张量 塑性应变张量 内变量 本步应变增量 这里上一步中的物理量用下标
n‑1标出，本步骤中的物理量在下标n中省略。塑性问题的数值计算常用返回映射算法，其由两个步骤组成，弹性试验步假设材料具有弹性，试应力 Cijkl表示弹性
张量， 表示上一步应力张量，Δεkl表示本步应变增量。此时的应力偏离该步骤的真实屈服面，然后使用径向返回算法进行塑性校正，将试应力拉回到屈服面上。对于J2流动理论的特殊情况，只有一个内部变量 屈服条件是
[0033]
[0034] 其中sij是偏应力， 是等效塑性应变，σy是当前屈服应力。径向返回算法一般使用循环迭代判断的方法实现，第k个迭代步的应力和塑性应变的增量为
[0035]
[0036]
[0037] 这里nij是屈服面的法向方向，且
[0038]
[0039] 其中，G是剪切模量，Ep是塑性模量，当|fk|＜TOL时迭代结束，其中TOL为给定的迭代收敛容差。
[0040] 根据Johnson和Cook提出的有关应变硬化、应变率和热软化的流动应力模型，则应变率效应必须包含在冲击和爆炸问题的本构模型中，
[0041]
[0042] 其中A、B、C、m1和m2是Johnson‑Cook模型的材料常数。 是塑性应变率， 是有效塑性应变率，Tr和Tm是室温和熔化温度，T是当前温度。Johnson‑Cook模型还定义了物质点随时间的累积损伤，
[0043]
[0044] 这里 为等效塑性应变增量，当D＝1时物质点发生断裂。 是断裂的等效应变，[0045]
[0046] 其中σm是三个法向应力的平均值，是Mises等效应力，D1、D2、D3、D4和D5是Johnson‑Cook损伤的材料常数。可以通过损伤来折减应力，当损伤参数D达到1时，发生断裂。断裂发生后，偏应力置零；如果球应力为正，则再将球应力置零。
[0047] 最后，导出时间间断态基近场动力学(SBPD‑TD)的基本格式。在时间间断伽辽金方法中，允许未知场函数在相邻时间间隔之间间断。函数w(t)在时间域的阶跃定义为，[0048]
[0049] 其中 是间断算子符号。对于任意的时间间隔In，其初时刻和末时刻分别为tn和tn+1，位移u和速度v可以分别基于三次Hermite函数和线性函数进行插值，
[0050]
[0051] 其中，vn和vn+1是本时间步tn和tn+1时刻的速度，un和un+1是本时间步tn和tn+1时刻的位移，且
[0052]
[0053] 时间步长Δt＝tn+1‑tn，I是单位矩阵。类似地，假设外力密度b在每个时间间隔内线性变化，内力密度f可以通过Hermite函数独立插值，
[0054]
[0055] 其中，bn和bn+1是本时间步tn和tn+1时刻的外力密度，fn和fn+1是本时间步tn和tn+1时刻的内力密度， 和 分别表示内力密度在本时间步tn和tn+1处的时间导数。
[0056] 将动量方程(7)写作
[0057]
[0058] 这里，
[0059]
[0060]
[0061] 其中Nxdim为Nx个物质点的自由度数。
[0062] 内力密度f可以线性化表示为
[0063] f＝Ku(23)虽然物质点的位移和速度是独立插值的，但它们必须满足如下约束条件，
[0064]
[0065] 运动方程，约束条件，连同时间间隔In内位移和速度的间断方程一起构成积分弱形式，
[0066]
[0067] 从δvn、δvn+1、δfn和δfn+1的独立变分得到SBPD‑TD的基本公式，即，[0068]
[0069]
[0070]
[0071]
[0072] 这里，ρ是密度， 是上一步计算出的末时刻位移， 是上一步计算出的末时刻速度， 是上一步计算出的末时刻内力，且
[0073]
[0074] 以上的求解格式是本发明提出的一种新的逐步积分法，用以替代求解动量方程(7)常用的中心差分法(CDM)，从而控制由中心差分法带来的应力解中的数值振荡。
[0075] 根据上述理论推导得出的时间间断态基近场动力学的基本格式，采用显式求解策略，再结合经典材料本构计算方法即可实现本发明提出的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，结合图1所示本发明提出的SBPD‑TD方法计算循环示意图，其具体实施步骤如下，
[0076] 步骤1：建立离散物质点模型，定义材料参数，根据离散的物质点位置坐标和预设的邻域大小判断每个物质点的邻域范围，建立物质点之间的相互作用关系；
[0077] 步骤2：定义边界条件和时间步长，根据位移场和速度场利用态基近场动力学方法计算结构的内力和内力的时间导数，引入经典的塑性本构及损伤演化法则进行非线性分析，，更新应力和塑性应变；
[0078] 步骤3：通过实施交错迭代求解策略，利用式(26)‑(29)进行显式时程积分，更新位移和速度，当迭代收敛时获得本步的位移解和速度解，返回步骤2进行下一时间步计算。
[0079] 在步骤1中，为了更加精确地贴合复杂结构的表面，可以使用有限单元法中的六面体网格进行离散，单元的形心作为物质点坐标，单元的体积作为物质点体积。作为一种非局部模型，近场动力学的物质点与其邻域内所有的物质点都具有相互作用。为了减少计算规模，可以基于胞元构建邻域。如图3(a)所示，物质点x的视界定义为以δ为特征的范围内所有物质点。如图所示划分胞元后，只需搜索x所在胞元及其周围8个胞元中的物质点即可得到物质点的邻域。(对于三维情况，需要搜索27个单元格，如图3(b)所示。)值得注意的是，由于接触邻域会随时间变化，这种快速邻域搜索算法可以大大减少接触模拟的时间。
[0080] 在步骤2中，为了抑制非常规态基近场动力学中固有的零能模式，本发明引入一种基于键的变形梯度修正方法，定义变形状态的修正，
[0081] zpq＝Ypq‑Fpξpq (31)将变形状态的非均匀部分添加到力状态中，
[0082]
[0083] 其中g是正常数，c0是键基近场动力学中的微模量，ω0定义为
[0084]
[0085] 在步骤2中，使用径向返回算法更新应力和塑性应变。J2流动理论的屈服面在π平面内是一个圆，其法向沿圆的径向。采用关联塑性流动法则，塑性流动方向为，
[0086]
[0087] 其中， 这里 表示试应力，
[0088] 在步骤2中，对于非线性分析，本构模型可以用率形式表示。柯西应力张量σij的物质时间导数与Jaumann应力率 有关，
[0089]
[0090] 其中Ωij是旋率张量，它是速度梯度L的反对称部分，
[0091]
[0092] 对于各向同性弹性模型，Jaumann应力率可以从变形张量率 导出，
[0093]
[0094]
[0095] 其中Cijkl是弹性张量。
[0096] 在步骤2中，可以引入两种损伤的形式。一种是基于键的损伤，一种是基于点的损伤。对于基于键的损伤，可以通过如下的三种方式判断键的断裂，
[0097] (1).定义物质点p和q之间键的伸长率为
[0098]
[0099] 其中ξpq是物质点p和q在参考构型中的相对位置，ηpq是物质点p和q之间的相对位移。键的断裂状态μ可由下式给出，
[0100]
[0101] 其中s0表示键的临界伸长率。
[0102] (2).键的断裂状态也可以直接通过应力或应变(如临界等效塑性应变 )来判断，[0103]
[0104] (3).键的渐进性延性损伤也可以定义为，
[0105]
[0106] 其中Dc≤1是渐进性延性损伤材料参数(Dc＝1对应于键的突然断裂)；Dξ是键的损伤参数，可以通过两端物质点的Johnson‑Cook损伤加权平均获得。
[0107] 基于上述三种方法，在基于键的损伤模型中，近场动力学的标量影响函数ωpq需要修改为ωpq＝μωpq。此外，还需引入了一个标量函数 来描述物质点xp的失效程度，[0108]
[0109] 对于基于点的损伤，可以借鉴经典的连续介质损伤模型。例如，通过Johnson‑Cook模型定义物质点的损伤后，可以通过损伤来降低应力。当损伤参数D达到1时，发生断裂。断裂发生后，偏应力s设置为0。如果计算的球应力为正，则将其重置为零。
[0110] 在步骤3中，式(26)‑(29)的显式时程积分方案使用不动点迭代的方法完成。其具体实施过程使用伪代码的形式给出：
[0111] (1)初始化变量，令k＝0。
[0112] (2)若k＜kmax，开始循环；
[0113] (2.1)根据近场动力学方法，使用速度 和位移un、 分别计算内力密度的时间导数 和
[0114] (2.2)根据近场动力学方法，使用位移 计算内力密度
[0115] (2.3)根据式(26)‑(29)，计算
[0116] (2.4)若 转进入步骤(3)，否则k＝k+1，返回步骤(2)继续循环。
[0117] (3)更新速度 位移 上述的实施步骤中n表示时间步。
[0118] 结合图1所示的SBPD‑TD的操作流程图，本发明所提出的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法(SBPD‑TD)的具体实现过程将通过如下的伪代码形式展示：
[0119] (1)建立离散的物质点模型，构建邻域，定义材料参数(密度ρ、弹性模量E、泊松比pν；塑性模量E、Johnson‑Cook参数(A、B、C、m1、m2)、初始屈服应力σy0、临界伸长率s0、临界等效塑性应变 渐进性延性损伤材料参数Dc等)，施加载荷并初始化速度位移等各变量。
[0120] (2)时间步循环n＝0,1,…,Nt。
[0121] (2.1)进行变量更新，令k＝0；
[0122] (2.2)若k＜kmax，开始循环；
[0123] (2.2.1)使用速度 和位移un、 分别计算内力密度的时间导数 和以及内力
[0124] (2.2.2)根据式(26)‑(29)，更新
[0125] (2.2.3)若 进入步骤(2.3)；否则k＝k+1，返回步骤(2.2)继续循环。
[0126] (2.3)更新速度 位移 若n＜Nt，则n＝n+1，返回步骤(2)；
[0127] 否则时间步循环结束，进入步骤(3)。
[0128] (3)输出计算文件，进行后处理。
[0129] 本发明的有益效果：
[0130] 采用本发明提供的技术方案与现有技术相比，具有如下显著效果：
[0131] (1)本发明提供的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，为结构冲击断裂分析提供了一种简便的数值仿真方法，并且拓宽了近场动力学方法在固体材料断裂领域的应用范围。通过采用近场动力学模型可以有效克服采用传统连续介质模型难以准确捕捉裂纹扩展路径的不足，可以较为简单地处理复杂裂纹扩展问题，例如裂纹交叉、分岔及裂纹在三维空间自由扩展等。此外，使用无网格方法离散的近场动力学模型还降低了数值实施过程的复杂度。
[0132] (2)本发明提供的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，将时间间断的思想有效地引入了近场动力学框架中，为结构冲击问题中的应力分析提供了一种高精度的数值仿真方法。通过使用时间间断方法，有效地控制了传统近场动力学求解时由于采用中心差分法(SBPD‑CDM)而带来的高频数值振荡，为精确预测损伤和裂纹的起始和演化奠定了坚实的基础。本发明所提出的方法作为一种新的时程积分方法，可以方便地替换传统时程积分方法，减少了程序编制的难度。
[0133] (3)本发明提供的结构冲击弹塑性断裂分析的时间间断态基近场动力学方法，尤其提供了一种有效的结构冲击破坏显式弹塑性分析的数值仿真方法，将经典材料本构模型有效地引入了态基近场动力学，克服了传统键基近场动力学的材料参数限制和只能分析弹脆性材料的缺点。通过引入多种损伤判据，使本发明可以有效预测拉伸型裂纹和剪切型裂纹的扩展路径，并准确识别裂纹类型的转化。采用的近场动力学法因其在处理大变形及极端变形时具有独特优势有效克服了基于传统有限单元法模拟有限变形弹塑性断裂破坏时因发生网格畸变而造成数值精度严重降低的不足。
附图说明
[0134] 图1为本发明的一种结构冲击弹塑性断裂分析的时间间断态基近场动力学方法(SBPD‑TD)的操作流程图；
[0135] 图2为本发明的非常规态基近场动力学模型的初始和当前构型示意图；
[0136] 图3为本发明的一种基于胞元的(a)二维和(b)三维快速邻域搜索算法示意图；
[0137] 图4为本发明的实施例1弹塑性杆中应力波的传播结构及边界条件示意图；
[0138] 图5为本发明的实施例1弹塑性杆在(a)0.055s和(b)0.135s时分别用SBPD‑TD和SBPD‑CDM求得的应力解与有限元解的对比；
[0139] 图6为本发明的实施例2分别由(a)SBPD‑CDM和(b)SBPD‑TD求得的应力解和断裂状态；
[0140] 图7为本发明的实施例2在(a)0.12s和(b)0.13s时分别由SBPD‑TD和SBPD‑CDM求得的应力解和断裂状态；
[0141] 图8为本发明的实施例3模型尺寸及边界条件示意图；
[0142] 图9为本发明的实施例3在40μs时由SBPD‑TD和SBPD‑CDM得到的Mises应力分布和断裂情况；
[0143] 图10为本发明的实施例3分别由SBPD‑TD和SBPD‑CDM计算得到的最终损伤形貌；
[0144] 图11为本发明的实施例4两变形体冲击碰撞模型及边界条件示意图；
[0145] 图12为本发明的实施例4的(a)铜制方板和(b)铝制圆柱在60μs时的Mises应力分布和碰撞形貌；
具体实施方式
[0146] 下面结合附图和实施例对本发明的性能做出进一步详细说明。以下实施例用于说明本发明，但不能用来限制本发明的适用范围。
[0147] 为了使本发明的目的、技术方案和具体实施效果展示更加清晰明了，下面通过四个具体实施例结合附图4～12对本发明提出的SBPD‑TD方法的准确性和有效性作进一步的详细说明。首先，通过两个一维标准算例分别验证本发明所提出方法在控制数值振荡和准确预测裂纹位置的准确性和有效性。然后，通过采用二维和三维算例说明本发明提出的SBPD‑TD方法抑制高维问题应力解中数值振荡的能力。
[0148] 实施例1：弹塑性杆中应力波的传播(附图4～5)
[0149] 本实施例为标准数值算例，考虑了一维结构中应力波的传播问题，其结构尺寸和边界条件如图4所示。杆的长度为10m，被离散成5000个物质点。杆的右端自由，其左端受到冲击荷载，
[0150]
[0151] 使用Johnson‑Cook本构模型，杨氏模量E＝10000Pa，泊松比ν＝0，质量密度ρ＝3
1kg/m 。Johnson‑Cook模型的材料常数为A＝5Pa，B＝600Pa，m1＝0.7，C＝0。图5说明了在ta＝0.055s和tb＝0.135s分别使用时间间断态基近场动力学SBPD‑TD和近场动力学耦合中心差分法SBPD‑CDM获得的应力波分布。由于塑性模量随应力状态而变化，应力波在杆中的传播变得相当复杂，很难得到解析解。因此，本实施例考虑与采用中心差分法时程积分方案的有限元法FEM‑CDM进行比较。数值结果表明，SBPD‑TD可以准确模拟应力波传播的高梯度特征并控制由中心差分法带来的虚假数值振荡。此实施例有效地说明了本发明提出的时间间断态基近场动力学方法在应力波传播模拟时的有效性和准确性，并验证了弹塑性程序的正确性。
[0152] 实施例2：弹塑性杆断裂位置的预测(附图6～7)
[0153] 本实施例为弹塑性杆受冲击荷载时断裂位置的预测算例。近场动力学作为一种非局部理论，可以在一定程度上克服传统连续介质力学损伤研究中遇到的奇异性难题。因此，近场动力学在处理断裂问题方面具有显著优势。下面通过一个数值算例来说明时间间断态基近场动力学在预测裂纹位置方面的优势和必要性。首先考虑一根可受压但只能承受较低拉应力的杆，模型的尺寸与边界条件如图4所示。杆的长度为10m，被离散成5000个物质点。杆的右端自由，其左端受到冲击荷载如式(44)所示。使用Johnson‑Cook本构模型，杨氏模量
3
E＝10000Pa，泊松比ν＝0，质量密度ρ＝1kg/m。Johnson‑Cook模型的材料常数为A＝5Pa，B＝600Pa，m1＝0.7，C＝0。临界失效应力设置为σf＝1Pa。图6显示了分别使用SBPD‑CDM和SBPD‑TD获得的0.09秒时沿杆的损伤和应力分布。在应力波在自由端反射之前，拉应力不会出现在杆中。然而，如图6(a)所示，CDM给出的应力解具有虚假的拉应力，这种虚假的拉应力又导致了虚假的断裂。这个例子强调了时间间断方法对应力波传播及断裂判断研究的重要性。
[0154] 接着提高临界失效应力为σf＝3.2Pa，相对可受压但只能承受较低拉应力的假设，这个假设主要考察杆在发生自由端反射时的断裂行为。当应力波在右端自由界面反射时，压应力会转化为拉应力，这就是层裂发生的原因。图7(a)给出了在0.12s时SBPD‑TD和SBPD‑CDM计算的应力分布。此时，由于应力波的叠加，拉应力刚刚出现在杆内。由于SBPD‑CDM解出现剧烈的数值振荡，拉应力超过3.2Pa，导致杆发生了断裂。但根据ABAQUS的参考解，可以推断此时的应力不超过3.2Pa，SBPD‑CDM解的断裂是由于虚假数值振荡导致应力提前满足断裂准则。真正的断裂位置应该如图7(b)所示，即断裂首先发生在8.9m处，这也被ABAQUS的应力分析所证实。
[0155] 上述两算例给出的数值结果与理论和经典数值分析结果吻合良好，很好地验证了本发明所提出的时间间断态基近场动力学方法在弹塑性冲击断裂分析中的准确性和有效性，可以有效地捕捉应力的高梯度特性并抑制虚假的数值振荡，从而准确地判断了断裂的发生位置和时间，为后续高维结构分析奠定了良好的基础。
[0156] 实施例3：二维弹塑性板的冲击断裂问题(附图8～10)
[0157] 本实施例为二维冲击断裂试验的验证算例，该试验由Zhou等人完成。图8为模型尺寸及边界条件(长度单位为mm)，将板均匀地划分为100×200个物质点。其弹性模量E＝3
192GPa，泊松比ν＝0.3，质量密度ρ＝7830kg/m 。采用线性硬化模型，初始屈服应力σy0＝p
2GPa，塑性模量E ＝200GPa。冲击速度v＝25m/s。根据试验观测，当冲击速度20m/s＜v＜
29m/s时，在裂纹的初始扩展阶段，断裂的主要机制是材料的剪切行为，剪切带从缺口尖端沿几乎平行于冲击方向的方向传播；在剪切裂纹扩展一定距离之后，裂纹转化为拉伸型并从剪切带的尖端沿与冲击方向成大约30°的方向延伸。从图9可知，SBPD‑TD方法有效地抑制了SBPD‑CDM带来的严重数值振荡。板的最终的损伤形貌如图10所示。
[0158] 实施例4：两变形体冲击碰撞行为模拟(附图11～12)
[0159] 最后一个实施例将考察时间间断态基近场动力学在实际案例中的有效性。如图11所示，一个铝圆柱以v0＝500m/s的速度撞击一个铜制方板。圆柱的材质为6061‑T6铝，被分3
为21万个物质点，其弹性模量E＝71GPa，泊松比ν＝0.28，质量密度ρ＝2704kg/m。Johnson‑Cook模型参数A＝324MPa，B＝114MPa，m1＝0.42，C＝0.002。方板的材料为黄铜，被平均分为
3
80万个物质点。杨氏模量E＝115GPa，泊松比ν＝0.31，质量密度ρ＝9095kg/m 。Johnson‑Cook模型参数A＝206MPa，B＝505MPa，m1＝0.42，C＝0.01。图12给出了SBPD‑TD在60μs时计算的Mises应力分布以及两个变形体的相对位置，可以看出SBPD‑TD有效地控制了应力解中的数值振荡。
[0160] 综上所述，我们首先通过实施弹塑性杆中应力波传播的研究验证了本发明所提出的时间间断态基近场动力学方法SBPD‑TD在应力波传播分析中的重要性和准确性，说明了应力波数值振荡的程度及控制方法。随后，考察弹塑性杆的断裂行为预测算例，通过两种不同的断裂判据考察了杆内两种不同的应力波传播现象引发的断裂，有效说明了SBPD‑TD在弹塑性断裂分析中的正确性和必要性。接着，通过对二维冲击裂纹扩展试验的仿真，有效说明了SBPD‑TD在二维弹塑性波传播分析中的准确性，并再现了试验中的冲击裂纹扩展路径和裂纹类型的转化。最后，一个实际的三维接触碰撞问题被考察，说明了SBPD‑TD在抑制波传播过程中数值振荡的有效性，同时为SBPD‑TD在大规模工程实际问题的分析奠定了坚实的基础。因此，本发明所提出的时间间断态基近场动力学方法SBPD‑TD是一种极具发展前景的冲击断裂破坏分析数值算法。
[0161] 本发明的实施例是为了示例和描述起见而给出的，而并不是无遗漏的或者将本发明限于所公开的形式。很多修改和变化对于本领域的普通技术人员而言是显而易见的。选择和描述实施例是为了更好的说明本发明的原理和实际应用，并且使本领域的普通技术人员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。
