CN109117508B - 一种用于模拟金属冷喷涂的无网格数值计算方法 - Google Patents

Abstract

本公开提供了一种用于模拟金属冷喷涂的无网格数值计算方法，包括：设置初始无网格计算模型参数；生成虚拟粒子，用以施加边界条件；进行邻近粒子搜索；计算每一粒子的核函数、采用核梯度修正法计算核函数梯度值以及采用双线性JC本构方程计算内力；计算动量及能量变化率，并采用变光滑长度法更新粒子的光滑长度；进行时间积分，并输出金属冷喷涂的模拟结果。本公开用于模拟金属冷喷涂的无网格数值计算方法具有高精度，适用于高应变率下的金属冷喷涂模拟，计算稳定性高。

Description

一种用于模拟金属冷喷涂的无网格数值计算方法

技术领域

本公开属于计算力学技术领域，更具体的涉及一种用于模拟金属冷喷涂的无网格数值计算方法。

背景技术

含多相、多介质的金属冷喷涂问题是高性能增材制造技术(3D打印)研究的热点问题，涉及运动界面、移动边界、大变形、大密度比等复杂特征。冷喷涂颗粒在与基板的超高速撞击作用下发生相变、破碎、结合甚至熔化，很难用基于网格的传统方法进行高效、可靠的数值模拟。传统的网格类方法通常基于动网格技术或嵌套网格技术来解决网格变形问题，然而当所研究问题的边界变形过大时采用这两种技术仍可能会产生网格畸变，此时通常利用网格重构来消除畸变的网格，然而网格重构必然会降低精度和增大计算量。为此人们提出了无网格方法，无网格方法不需要网格，因此不会遇到网格畸变的问题，能够自然追踪运动界面，在处理大变形及流固耦合问题时有着特殊优势。但是现有的无网格数值模型都不能很好的模拟金属冷喷涂问题，因此有必要提出适用于模拟金属冷喷涂的无网格数值计算模型。

2010年Li等人发展了一种用于模拟金属冷喷涂的SPH(smooth particlehydrodynamics，光滑粒子流体动力学)计算方法。其存在以下缺陷：

(1)该方法精度偏低，只有零阶精度；

(2)该方法使用Johnson-Cook本构模型，不适用于高应变率下的金属冷喷涂模拟；

(3)该方法采用固定光滑长度，在模拟大变形问题时稳定性差，同时加剧了数值模拟误差；

总之，该方法无法有效模拟多颗粒下的金属冷喷涂问题。

发明内容

(一)要解决的技术问题

基于以上问题，本公开的目的在于提出一种用于模拟金属冷喷涂的无网格数值计算方法，在SPH方法的基础上进行改进，用于解决以上技术问题的至少之一。

(二)技术方案

为了达到上述目的，作为本公开的一个方面，提供了一种用于模拟金属冷喷涂的无网格数值计算方法，包括：

设置初始无网格计算模型参数；

生成虚拟粒子，用以施加边界条件；

进行邻近粒子搜索；

计算每一粒子的核函数、核函数梯度值以及内力；

计算动量及能量变化率，并采用变光滑长度法更新粒子的光滑长度；

进行时间积分，并输出金属冷喷涂的模拟结果。

在一些实施例中，根据粒子的密度变化率更新粒子的光滑长度，具体关系式如下：

式中，h_i和ρ_i分别为粒子i的光滑长度及密度。

在一些实施例中，在计算每一粒子的核函数、核函数梯度值之后还包括：采用核梯度修正法对核函数的梯度进行修正，计算过程如下：

式中W_ij为核函数，

为原有核函数梯度，

为修正后的核函数梯度，粒子i位置信息为(x_i，y_i，z_i)，粒子j位置信息为(x_j，y_j，z_j)，公式中x_ji＝x_j-x_i，y_ji＝y_j-y_i，z_ji＝z_j-z_i，V_j为粒子j的体积。

在一些实施例中，采用双线性Johnson-Cook本构模型进行内力计算，计算过程如下：

式中，T^*为温度相关无量纲量；

为有效塑性应变，

为有效塑性应变率，

为准静态参考应变率；A，B，n，m均为材料常数；C和C′分别为低于和高于临界应变率

下的应变率敏感性指数。

在一些实施例中，所述初始无网格计算模型参数包括实体粒子信息，该实体粒子信息包括冷喷涂颗粒、基板粒子的位置信息及物理信息。

在一些实施例中，在所述生成虚拟粒子，并施加边界条件的步骤中，生成计算所需要的虚拟粒子，布置固壁以及镜像边界虚拟粒子信息，包括：

生成Monaghan型虚拟粒子，为第1层边界粒子，等间距的布置在边界上，其对接近的实际粒子产生排斥力，施加固定壁面边界条件；以及

生成Libersky型虚拟粒子，为第2～3层边界粒子，由靠近边界位置的实体粒子关于边界线镜像得到，其和实体粒子速度方向相反，施加防止粒子穿透边界的镜像边界条件。

在一些实施例中，根据所述实体粒子信息及虚拟粒子信息，采用三维树形搜索法进行邻近粒子搜索，确定每一粒子的相互作用粒子。

在一些实施例中，根据连续性方程、动量方程及能量方程计算每一时间步的动量及能量变化率。

在一些实施例中，所述核函数为三次样条核函数，其表达式如下：

式中，R是基于光滑长度h进行无量纲化后的粒子间距，R＝r/h，r为粒子间实际距离；α_d是根据归一性条件确定的常数，其在一维、二维和三维中分别取为1/h，15/7h²和3/2h³。

在一些实施例中，采用蛙跳法或预估校正法进行时间积分。

(三)有益效果

(1)本公开用于模拟金属冷喷涂的无网格数值计算方法具有一阶精度，与现有模拟金属冷喷涂的SPH方法相比，本公开精度更高。

(2)本公开用于模拟金属冷喷涂的无网格数值计算方法可以准确描述金属材料在高应变率下的本构关系，适用于高应变率下的金属冷喷涂模拟。

(3)本公开用于模拟金属冷喷涂的无网格数值计算方法适用于对粒子高度聚集及分散的模拟，计算稳定性高。

附图说明

图1为本公开用于模拟金属冷喷涂的无网格计算模型流程图；

图2为本公开边界虚粒子施加示意图；

图3为本公开改进Johnson-Cook模型计算的屈服应力随应变率变化曲线及相关实验数据；

图4为本公开实施例1采用的单颗粒金属冷喷涂粒子布置图；

图5为本公开实施例1生成的虚粒子布置图；

图6为本公开实施例1后冷喷涂颗粒最终变形图；

图7为本公开实施例1某时刻系统等效应力云图；

图8为本公开实施例1某时刻颗粒等效应力云图；

图9为本公开实施例1某时刻系统温度云图；

图10为本公开实施例2采用的多颗粒金属冷喷涂粒子布置图；

图11为本公开实施例2生成的虚粒子布置图；

图12为本公开实施例2单层颗粒冷喷涂某时刻等效应力云图；

图13为本公开实施例2单层颗粒冷喷涂某时刻温度云图；

图14为本公开实施例2双层颗粒冷喷涂时刻1等效应力云图；

图15为本公开实施例2双层颗粒冷喷涂时刻1温度云图；

图16为本公开实施例2双层颗粒冷喷涂时刻2等效应力云图；

图17为本公开实施例2双层颗粒冷喷涂时刻2温度云图；

图18是现有技术中采用SPH方法模拟的多颗粒金属冷喷涂形态图。

<符号说明>

1-实体粒子；1′-基板粒子；2-Monaghan型虚拟粒子；3-Libersky型虚拟粒子。

具体实施方式

为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开作进一步的详细说明。

本公开提供了一种用于模拟金属冷喷涂的无网格数值计算方法，如图1所示，所述用于模拟金属冷喷涂的无网格数值计算方法包括以下步骤：

步骤1、首先设置初始无网格计算模型参数，包括冷喷涂颗粒及基板粒子的初始位置(x，y，z)信息和物理信息，所包含物理信息如速度v、密度ρ、质量m等。

步骤2、生成计算所需要的虚拟粒子，布置固壁以及镜像边界粒子信息(v、ρ、m等)，用以模拟边界条件以及防止粒子穿透。

如图2所示，第1层边界粒子为Monaghan型虚拟粒子，等间距的布置在边界上，该Monaghan型虚拟粒子对接近的实际粒子产生排斥力，施加固定壁面边界条件。同时边界外布置2～3层Libersky型虚拟粒子，该Libersky型虚拟粒子由靠近边界位置的实体粒子关于边界线镜像得到，有和实体粒子相同大小的物理量，速度方向相反，施加防止粒子穿透边界的镜像边界条件。

步骤3、根据步骤1、2中布置的粒子位置信息(x，y，z)及粒子物理信息(v、ρ、m等)，采用三维树形搜索法进行邻近粒子搜索，确定每一粒子(实体、虚拟粒子)的相互作用粒子，用以计算粒子(实体、虚拟粒子)在每一时间步的物理量。

步骤4、先计算每一粒子(实体、虚拟粒子)的核函数及其导数值，再采用核梯度修正法(KGC)对核函数的梯度进行修正，由此能够提高模拟的精度。

其中，核函数有很多种，这里仅示例性的给出三次样条核函数，但本公开核函数并不限于三次样条核函数。公式(1)是三次样条核函数的表达式。公式(1)中R是基于光滑长度h进行无量纲化后的粒子间距，R＝r/h，r为粒子间实际距离。α_d是根据归一性条件确定的常数，其在一维、二维和三维中分别取为1/h，15/7h²和3/2h³。

计算出核函数及其导数值之后，再采用核梯度修正法对核函数的梯度值进行修正。基于对SPH核近似的泰勒展开，可以得到公式(2)，其中O(h²)为泰勒展开的高阶余项。并将公式(2)与传统SPH求解场变量f(x)导数的公式(3)相结合，可以得到求解场变量导数的展开公式(4)。使用SPH粒子近似法得到求解场变量导数的粒子表达形式(5)，式中V_i，V_j分别为i，j粒子的体积。通过保证

和

使得对场变量导数的求解具有一阶精度。为此，对核函数的梯度

乘以一个局部可逆矩阵L(r_i)，如公式(6)和(7)所示。公式中x_ji＝x_j-x_i，y_ji＝y_j-y_i，z_ji＝z_j-z_i。将得到的新的核梯度

用于对场变量导数的计算中，可以大幅度提高计算精度。

步骤5、采用改进的双线性Johnson-Cook本构模型进行内力计算，此模型更适合高应变率下的内力计算，计算过程如下；

在偏应力的求解过程中，引入不受刚体转动影响的客观张量焦曼应力率

如公式(8)，式中S和

分别为偏应力及偏应力率，α，β，γ为张量指标。当一点的应力状态处于弹性阶段时，本构关系由公式(9)表示，式中

和

分别为速度梯度张量

的对称分量和反对称分量，如公式(10)所示。

为体积应变率。记弹性偏应力为S_e ^αβ，弹性偏应力率为

由式(8)和(9)可得公式(11)，其中G为剪切模量，δ为单位张量。弹性偏应力第二不变量J由公式(12)所示，材料的屈服应力为σ_Y。当J≤σ_Y ²/3和J＞σ_Y ²/3时，此时一点的应力状态分别处于弹性和塑性阶段，相应的偏应力可由公式(13)计算。

在模拟金属冷喷涂过程中，金属材料屈服应力主要受到应变率、应力硬化和热效应的影响。冲击动力学中常用的Johnson-Cook本构模型并不适合超高应变率的情况。本公开采用改进的双线性Johnson-Cook模型，如公式(14)所示。

其中

为温度相关无量纲量，T_room和T_melt分别为室温和金属熔化温度。

为有效塑性应变，

为有效塑性应变率，

为准静态参考应变率。A，B，n，m均为材料常数，可由实验获得。C和C′分别为低于和高于临界应变率

下的应变率敏感性指数，可由相应实验数据得到。

温度T和内能e之间的关系由公式(15)表示，式中C_V为定容比热容，下标“0”表示初始时刻物理量。计算中采用式(15)的一次时间导数形式进行温度更新，如公式(16)所示。其中de为内能变化量，上标“(n)”和“(n+1)”分别表示SPH计算中第n步和第n+1步物理量。当一点的应力状态进入塑性阶段时，需要单独计算偏应力做功对能量变化产生的影响即塑性功增量ΔW_p。第n步的塑性功增量

由公式(17)，(18)和(19)计算，式中

为中间步的密度，σ_eff和ε_eff分别为有效塑性应力和有效塑性应变。

e-e₀＝mC_V(T-T₀) (15)

步骤6、计算动量及能量变化率，用以求解下一步粒子(实体、虚拟粒子)的物理量，采用变光滑长度法更新粒子(实体、虚拟粒子)光滑长度，此方法能提高模拟的稳定性。具体更新方法如下；

在模拟金属冷喷涂过程中，超高速冲击下颗粒及基板粒子会出现聚集或是分散现象。因此需要在模拟过程中更新粒子的光滑长度，以保证每个粒子的周围既有足够的粒子又不会有过多的粒子与其相互作用。将粒子的光滑长度根据其密度变化率进行调整，调整如公式(20)所示。式中h_i和ρ_i分别为粒子i的光滑长度及密度，t为时间。

步骤7、采用蛙跳法或预估校正法等积分格式对所有物理量进行时间积分，更新粒子(实体、虚拟粒子)信息，一定时间步后输出最终模拟结果。

本公开用于模拟金属冷喷涂的无网格数值计算方法具有一阶精度，与现有模拟金属冷喷涂的SPH方法相比，本公开精度更高。具体说明如下：

若任意函数f(x，y，z)在点i(x，y，z)的某一邻域内具有n阶连续性，那么就可以在该点对函数f(x，y，z)进行n阶泰勒展开。为了方便，这里用f_i表示，f_i，α分别表示f(x，y，z)在i点处的场函数值及其在α方向的导数。保留函数f(x，y，z)在点i进行泰勒展开的一阶项，余项用r((x-x_i)²)表示，则可得下式：

公式(21)因保留了一阶项而具有一阶精度。在公式(9)两边同时乘以核函数及其在三个方向的导数，并在问题域内进行积分，可以得到下列方程组：

方程组(22)中Ω为相应的问题区域，为简便起见在下面公式中不再标出。将方程组(22)写为矩阵方程形式为：

则函数f及其导数可由下式求解：

W_i，x分别表示核函数在i点处对x方向的导数，dV表示相应的体积分。因为核函数的归一性及对称性，在内部求解区域∫W_idV＝1，∫W_i，xdV＝0，故可以忽略矩阵L^-1的第一行和第一列元素。进而将公式(24)和(25)写成相应的粒子求和形式可得：

上式中N为粒子i的邻近粒子数。比较公式(27)与公式(7)可知，L-¹和L(r_i)为同一个矩阵。公式(26)右端第一项可以写成

因此对核梯度

乘以一个可逆矩阵即公式(6)，可实现同公式(26)相同的计算结果。

在推导过程中只有公式(21)是近似表达式，其它公式的推导过程没有采用任何近似，因此公式(26)与公式(21)拥有相同的精度，根据公式(21)的精度可知本方法具有一阶精度。

本公开用于模拟金属冷喷涂的无网格数值计算方法可以准确描述金属材料在高应变率下的本构关系，采用改进的双线性Johnson-Cook本构关系适合计算材料高应变率下的内力。

传统的Johnson-Cook模型如公式(28)所示。该模型在从低至高的应变率范围内采用统一的公式计算屈服应力。图3为改进的双线性本构关系应力-应变率曲线以及相应的实验数据。在实际的数值模拟中，传统的Johnson-Cook模型在应变率小于10⁴s^-1时比较准确，而在金属冷喷涂过程中，材料的应变率可达到10⁷s^-1。这时Johnson-Cook模型将不再适用于模拟金属冷喷涂问题。而如图3所示，本公开采用改进Johnson-Cook模型由两阶段组成，在从低至高的应变率范围内和实验结果均比较吻合，因此其能适用于求解高应变率下的金属冷喷涂问题。

此外，本公开方法适用于对粒子高度聚集及分散的模拟，计算稳定性高。

下面结合实施例对本公开方法作详细说明。

实施例1

本实施例以单颗粒金属冷喷涂为例来说明本公开用于模拟金属冷喷涂的无网格数值计算方法。

步骤1

以单颗粒金属冷喷涂为例，基板的长L＝90μm，宽W＝90μm，厚T＝43μm，材料为铜。冷喷涂颗粒采用直径为16.5μm的球型铜颗粒，温度为350℃，以650m/s的初始速度冲击基板。选定初始时刻颗粒粒子间距dx₁＝dy₁＝dz₁及光滑长度h₁＝λdx₁。颗粒及基板粒子的密度均采用铜的密度ρ_cu＝8960kg/m³，粒子的质量m_cu＝ρ_cudxdydz。初始时刻基板粒子静止，冷喷涂颗粒粒子有着z方向的初始速度。图4给出了该模型的初始粒子布置图。

步骤2

在下边界处施加一层固体壁面边界粒子，边界外与实体粒子关于边界对称处生成两层镜像边界粒子。生成的虚拟粒子如图5所示。

步骤3

采用三维的树形搜索法进行邻近粒子搜索。通过粒子位置构造有序树，将问题域递归分割成一个个卦限，直到每一卦限只含有一个粒子。根据卦限分布及树型结构确定每一粒子的相互作用粒子。

步骤4

本算例采用三次样条核函数，具体表达式如公式(1)所示。计算每一粒子的核函数及其导数值，同时采用核梯度修正法(KGC)即公式(7)计算每一粒子新的核梯度值，并用新的核梯度替代原有核梯度进行后续计算。

步骤5

采用公式(13)计算离散SPH粒子在每一时间步内所受的偏应力，式中的屈服应力由双线性Johnson-Cook本构方程计算，即公式(14)。在一点的应力状态进入塑性阶段后，单独计算偏应力做功对能量变化产生的影响即塑性功增量

由此完成每一时间步的内力计算。

步骤6

根据连续性方程、动量方程及能量方程计算每一时间步的动量及能量变化率，用以求解下一步粒子的物理量。计算完成后根据公式(20)更新每个粒子的光滑长度。

步骤7

采用蛙跳法进行时间积分，更新粒子信息，每固定时间步数后输出模拟结果，包括粒子的位置、速度、密度、压力、温度、等效应力以及系统能量等信息。在冷喷涂颗粒达到稳定状态即不再运动后，输出最终模拟结果。图6为冷喷涂颗粒最终变形图，可以看到本公开的的数值模拟结果同实验结果非常接近。图7、图8分别为模拟过程中某时刻系统及颗粒的等效应力云图。图9为该模拟时刻下系统及颗粒的温度云图，图标为摄氏温度。由图可见，本模型能得到光滑的应力分布，在颗粒与基板界面处有少部分粒子超过熔点，而大部分界面粒子及内部粒子的温度低于熔点。这是理想的金属冷喷涂结果，即有大的塑性变形但没有熔化，因此不会产生由金属熔化、凝固引起的缩孔、偏析、裂纹等重要冶金问题。由此可见本公开的模型可准确模拟金属冷喷过程。

实施例2

本实施例以多颗粒金属冷喷涂说明本公开用于模拟金属冷喷涂的无网格数值计算方法。

步骤1

以单层及双层多颗粒金属冷喷涂为例。在单层颗粒冷喷涂中，铜基板的边长L＝110μm，宽W＝110μm，厚T＝43μm，基板上方布置6个冷喷涂颗粒。双层颗粒冷喷中，铜基板的边长L＝100μm，宽W＝100μm，厚T＝50μm，基板上方布置两层共8个冷喷涂颗粒。冷喷涂颗粒采用球型铜颗粒，以高速冲击基板。初始时刻模拟粒子间距为dx₁＝dy₁dz₁以及光滑长度h₁＝1.5dx₁。颗粒及基板粒子的密度均采用铜的密度ρ_cu＝8960kg/m³，粒子的质量m_cu＝ρ_cudxdydz。初始时刻基板粒子静止，冷喷涂颗粒粒子有着z方向的初始速度。图10给出了该模型单层颗粒(左图)及双层颗粒(右图)冷喷涂的初始粒子布置图。

步骤2

如图11所示，该问题在下边界处施加一层Monaghan型边界虚粒子，边界外与实体粒子关于边界对称处生成两层Libersky型边界虚粒子。

步骤3

采用三维的树形搜索法进行邻近粒子搜索，确定每一粒子的相互作用粒子。

步骤4

本算例采用三次样条核函数，具体表达式如公式(1)所示。采用核梯度修正法(KGC)即公式(7)计算每一粒子新的核函数导数值，用其替代原有核梯度进行后续计算。

步骤5

根据双线性Johnson-Cook本构方程即公式(14)计算SPH粒子在每一时间步内所受的偏应力。在一点的应力状态进入塑性阶段后，单独计算塑性功增量

完成内力计算。

步骤6

计算每一时间步的动量及能量变化率，计算完成后根据公式(20)更新每个粒子的光滑长度。

步骤7

采用蛙跳法进行时间积分，每固定时间步数后输出模拟结果，包括粒子的位置、速度、密度、压力、温度、等效应力以及系统能量等信息。图12、图13分别为单层颗粒金属冷喷过程中某时刻系统和颗粒的等效应力及温度分布云图，温度云图图标为摄氏温度。图14、图15、图16、图17分别为双层颗粒金属冷喷过程中两个时刻下系统和颗粒的等效应力及温度分布云图。

图18为现有技术中模拟的多颗粒金属冷喷涂结果，不同颜色代表不同冷喷涂颗粒。由图可见，在现有技术的模拟中，冷喷后金属颗粒内部散开，没有很好的成型，同时金属颗粒之间没有明显的结合。相比之下，本模拟能得到光滑的应力分布，冷喷涂颗粒发生大的塑性变形，实现不同颗粒冶金结合在一起。同时大部分界面模拟粒子及内部粒子的温度低于熔点，因此不会产生由金属熔化、凝固引起的缩孔、偏析、裂纹等重要冶金问题，这是理想的金属冷喷涂结果。由此可见，本公开可准确模拟金属冷喷涂过程。

在此处所提供的说明书中，说明了大量具体细节。然而，能够理解，本公开的实施例可以在没有这些具体细节的情况下实践。在一些实例中，并未详细示出公知的方法、结构和技术，以便不模糊对本说明书的理解。

类似地，应当理解，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开的示例性实施例的描述中，本公开的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开解释成反映如下意图：即所要求保护的本公开要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，公开方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本公开的单独实施例。

还需要说明的是，本文可提供包含特定值的参数的示范，但这些参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应值。实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向，并非用来限制本实用新型的保护范围。此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。

以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims (6)

1.一种用于模拟金属冷喷涂的无网格数值计算方法，包括：

设置初始无网格计算模型参数；

生成虚拟粒子，用以施加边界条件；

进行邻近粒子搜索；

计算每一粒子的核函数、核函数梯度值以及内力；

计算动量及能量变化率，并采用变光滑长度法更新粒子的光滑长度；

进行时间积分，并输出金属冷喷涂的模拟结果；

其中，根据粒子的密度变化率更新粒子的光滑长度，具体关系式如下：

式中，h_i和ρ_i分别为粒子i的光滑长度及密度，t为时间；

在计算每一粒子的核函数、核函数梯度值之后还包括：采用核梯度修正法对核函数的梯度进行修正，计算过程如下：

式中W_ij为核函数，

为原有核函数梯度，

为修正后的核函数梯度，粒子i位置信息为(x_i，y_i，z_i)，粒子j位置信息为(x_j，y_j，z_j)，公式中x_ji＝x_j-x_i，y_ji＝y_j-y_i，z_ji＝z_j-z_i，V_j为粒子j的体积；

采用双线性Johnson-Cook本构模型进行内力计算，计算过程如下：

式中，T^*为温度相关无量纲量；

为有效塑性应变，

为有效塑性应变率，

为准静态参考应变率；A，B，n，m均为材料常数；C和C′分别为低于和高于临界应变率

下的应变率敏感性指数。

2.根据权利要求1所述的用于模拟金属冷喷涂的无网格数值计算方法，其中，所述初始无网格计算模型参数包括实体粒子信息，该实体粒子信息包括冷喷涂颗粒、基板粒子的位置信息及物理信息。

3.根据权利要求2所述的用于模拟金属冷喷涂的无网格数值计算方法，其中，在所述生成虚拟粒子，用以施加边界条件的步骤中，生成计算所需要的虚拟粒子，布置固定壁面以及镜像边界虚拟粒子信息，包括：

生成Monaghan型虚拟粒子，为第1层边界粒子，等间距的布置在边界上，其对接近固定壁面的实际粒子产生排斥力，施加固定壁面边界条件；以及

生成Libersky型虚拟粒子，为第2～3层边界粒子，由靠近边界位置的实体粒子关于边界对称生成，其和实体粒子速度方向相反，施加防止粒子穿透边界的镜像边界条件。

4.根据权利要求3所述的用于模拟金属冷喷涂的无网格数值计算方法，其中，根据所述实体粒子信息及虚拟粒子信息，采用三维树形搜索法进行邻近粒子搜索，确定每一粒子的相互作用粒子。

5.根据权利要求1所述的用于模拟金属冷喷涂的无网格数值计算方法，其中，根据连续性方程、动量方程及能量方程计算每一时间步的动量及能量变化率。

6.根据权利要求1所述的用于模拟金属冷喷涂的无网格数值计算方法，其中，采用蛙跳法或预估校正法进行时间积分。

Images