本发明属于计算力学与数值仿真，具体涉及一种基于无网格法插值技术的冲击大变形过程模拟方法及应用。
背景技术：
1、在过去的几十年中，传统网格类方法飞速发展，已经被广泛地应用于科学和工程问题的研究中。但是受限于固有的网格属性，这类方法在解决冲击大变形问题的时候往往会遇到困难，比如拉格朗日网格（如有限元法）受网格畸变困扰，无法很好地解决超大变形、裂纹扩展等问题，欧拉网格（如有限差分法、有限体积法）很难追踪易变形边界和运动物质交界面。
2、无网格法不需要进行网格离散，在一定程度上解决了因网格依赖带来的部分计算难题，有利求解冲击大变形过程，在侵爆问题上应用广泛。
3、目前常用的无网格方法有光滑粒子流体动力学法（smoothed particlehydrodynamics, sph），无单元伽辽金法（element-free garlerkin, efg），再生核粒子法（reproducing kernel particle method, rkpm）以及物质点法（material point method,mpm）等。大量研究者使用以上方法及其改进方法对冲击大变形问题进行了模拟，虽取得了丰富研究成果，但仍然无法解决准确设置位移边界条件、计算效率不高、拉应力不稳定及有效计算带摩擦的动态接触等问题，同时也缺乏严密的收敛性与误差理论分析。最优运输无网格方法（optimal transportation meshfree, otm）是针对动态冲击、侵彻、爆炸等问题提出的显式增量更新拉格朗日无网格方法。otm方法在处理冲击大变形问题时具有以下几个优势：
4、1.该方法采用局部最大熵插值函数（local maximum entropy， lme），克服了一般无网格法中插值函数（如sph中的核函数近似和efg中的原始移动最小二乘插值）不满足kronecker-delta属性的本质缺陷，解决了传统无网格法难以准确施加位移边界条件的问题；
5、2.该方法采用节点和物质点组合离散的方式，以物质点作为积分点有效避免了计算结果在拉伸载荷下的不稳定性；
6、3.该方法采用基于能量准则的裂纹扩展算法，以物质点的材料能量释放率作为失效判据，能够很好地模拟材料的损伤情况，克服了传统裂纹扩展数值方法裂纹扩展路径网格相关、不收敛、无法清晰表征材料真实变形和失效物理机制的缺点。
7、与传统网格类方法相比，无网格方法的一个主要区别在于插值函数的处理。传统网格方法中的插值函数与单元捆绑，以具体的函数形式表征，可以通过函数赋值直接求得；而无网格方法的插值函数，在每次更新邻域粒子后，需通过拟合或优化重新计算，计算开销较大。因此，在无网格方法的计算流程中，高效稳定的插值技术显得尤为重要。otm方法使用的局部最大熵插值函数需通过求解凸优化问题来获得，但使用一般算法求解在模拟大变形问题时效率低且不稳定。例如，当材料变形过大时，邻域粒子减少，单一时间步内往往无法在物质点上求得收敛的形函数值，导致出现非物理的计算结果，进一步引发计算程序的崩溃。因此，本专利致力于完善局部最大熵插值函数的求解过程，解决otm方法中插值函数的瓶颈问题，提供高效稳定的仿真策略，以应对冲击大变形问题。
8、局部最大熵插值函数 n(x)（在凸优化问题中以概率 p(x)表征）可以通过求解一个凸优化问题获得：
9、，
10、其中x是待求插值函数的物质点坐标，应位于由点集 x构成的凸集conv x的内部， x a( a=1,…, n)是节点坐标，p(x)=( p1(x) , p2(x) ,… p n(x)) t，是待求的插值函数， β是帕雷托最优化参数，可以控制形函数的局部性，是一个给定的常数。空间中任意变量 u h(x)可以用节点值 u a= u a(x a)插值得到：
11、 ，
12、使用凸优化理论进行严密的推导，凸优化问题的对偶问题等价于：
13、，
14、代表以 λ为变量最小化函数的值，即为待求问题的目标函数；即为最终所需求解的无约束凸优化问题，此时插值函数p(x)可以表示为：
15、 ，
16、其中，，为统计力学意义上的配分函数，为范数运算符，此处一般使用2-范数。可以证明，目标函数 f( λ)的黑塞矩阵j( λ)是正定的，因此 f( λ)是r d上的严格凸函数，其最小值存在且唯一。
17、由于目标函数的形式过于复杂，无法直接通过极值定理来求解析解。一般地，对于这种无约束凸优化问题，可以通过牛顿法来迭代求解 f( λ)的最小参数值。牛顿法是一种基于二阶泰勒展开求解函数极值点的下降类方法，该类方法在每个计算时间都可以找到一个趋于极值的下降方向，依循该方向，总是可以得到一个比当前计算点更接近极值的点。以第k个计算时间步为例，选择不同的搜索算法可以得到不同的时间步长 t k，结合下降方向可以算出步进的计算点，其中，r和j分别是目标函数 f( λ)的梯度和黑塞矩阵。
18、但是，牛顿法的局部收敛性不可避免地带来了两个问题：
19、第一，在目标函数非凸的情况下，牛顿法常常失效，体现在空间中某些点的黑塞矩阵j的行列式为0，无法求逆，故无法求牛顿方向，因此在凸优化理论中往往会检验目标函数的黑塞矩阵是否正定；
20、第二，牛顿法的收敛性依赖初值的选择，在远离极值点的区域，牛顿法无法达到二阶收敛速度，不合适的初值会引发解的振荡甚至发散，因此在应用中常常会通过lipschitz连续条件来检验目标函数局部点的二次性态。
21、所要求解的公式就会经常性地遇到以下问题：
22、一个具体的局部最大熵形函数的求解问题：如附图1所示，二维点集 x取为{(0,0),(1,0), (1,1) , (0,1) }，计算点x取为(0.8,0.2)，参数取为1.6，可以画出的二维曲面（附图2），不难看出，该图像确实是r2上的严格凸函数。进一步地，该曲面在临近极值点处是光滑的弧面，呈现明显的二次性，而远离极值点的区域，被划分为了四块明显的平面，即当|| λ||→∞时， f( λ)趋近于线性函数，此时|| λ||=0，即黑塞矩阵是奇异的，行列式为0（附图3），不存在逆矩阵，将无法直接求得牛顿方向。更明确地说，当把初始值选在这些区域内，或迭代值落入这些区域内时，传统牛顿法将失去作用。
23、为解决上述问题，人们提出不同的办法来修正黑塞矩阵，其中比较著名的有阻尼牛顿法和拟牛顿法。阻尼牛顿法引入了阻尼参数，亦称惩罚因子，将原黑塞矩阵j修正为j+ μi，保证了下降方向的存在。当 μ较大时，该方法实际上是将牛顿法退化为了只有一阶收敛速度的最速下降法，全域使用同一个常数 μ，极大地影响着计算效率。拟牛顿法的基本思想是“对黑塞矩阵（或逆）做近似”，使其比最速下降法快，又比牛顿法计算简单，且整体收敛性好。由拟牛顿条件可以引申出不同的拟牛顿法，其中使用得比较广泛的有dfp算法和bfgs算法，但计算效率依旧存在问题。
技术实现思路
1、本发明的目的在于提供一种基于无网格法插值技术的冲击大变形过程模拟方法，以解决局部最大熵形函数的求解冲击问题的过程中，通过牛顿法来迭代求解f(λ）的最小参数值时可能会遇到失效、计算效率低等问题。
2、为了实现上述目的，本发明的技术方案如下：
3、本发明涉及一种基于无网格法插值技术的冲击大变形过程模拟方法，包括以下步骤：
4、s1.设置冲击大变形问题域，设置全局参数和包括材料模型在内的模型参数；
5、s2.将问题域离散为一组物质点集和一组节点集，为每个物质点初始化邻域节点集，为节点初始化物理变量，物理变量包括节点集中质量、线性动量和节点力；
6、s3.离散时间域，以时间增量[ t n, t n+1]为单位循环进行s4-s8；
7、s4.对于每个物质点，在邻域节点集上构建目标函数 f( λ)，并采用正则化牛顿法对问题域中的变量 λ进行迭代求解，保存最终变量 λ*，计算局部最大熵插值函数及其导数；
8、s5.基于局部最大熵插值函数更新节点集中质量，更新下一时刻的节点坐标，基于下一个节点的节点坐标以及最大熵插值函数的导数计算局部变形梯度增量，基于局部变形梯度增量计算下一时刻的局部变形梯度；
9、s6.将s5中求得的局部变形梯度增量和局部变形梯度代入材料模型中，求得物质点的应力；
10、s7.综合s6中求得的物质点的应力和边界条件，更新节点力和线性动量，判断是否计算至最后一个时间步，若不是，进入s8，若是，进入s9；
11、s8.更新物质点坐标和邻域节点集，返回s4；
12、s9.完成冲击过程模拟，输出结果。
13、优选地，所述s4中采用正则化牛顿法对问题域中的变量 λ进行迭代求解的具体步骤为：
14、s4.1.根据初始化阶段确定的节点集和待求物理变量的初始值 λ0，设定精度阈值 δ，令待求物理变量 λ=λ0；
15、s4.2.计算每个节点上的配分函数值 z a，计算公式为：
16、 （1），
17、公式中， β为帕雷托最优化参数，x表示计算的物质点，x a表示节点集{x a}中的任意一个节点，||·||表示范数运算；
18、对各节点的配分函数值进行求和得到；
19、s4.3.计算形函数 n a，计算公式为：
20、 （2） ；
21、s4.4.对目标函数 f( λ)进行正则化修正，修正后的目标函数为，其表达式为：
22、（3），
23、其中， λ为正则化修正点，即当前时刻的计算点， v为正则化修正函数的自变量，▽ f( λ)表示目标函数的梯度向量；
24、s4.5.计算目标函数的梯度向量▽ f( λ)，计算公式为：
25、（4），
26、公式中， n表示节点的数量， a表示节点索引号；
27、s4.6.对梯度向量进行范数运算获得精度，并与设定的精度阈值 δ进行比较，若计算所得的进度不大于精度阈值，输出此时的变量作为最终变量 λ*，进入s4.9，否则进行s4.7- s4.8；
28、s4.7.计算目标函数的黑塞矩阵h（ λ），基于黑塞矩阵及梯度向量计算正则化牛顿下降方向，
29、黑塞矩阵的计算公式为：
30、（5），
31、公式中，r(x, λ)表示梯度向量，i表示单位矩阵；
32、正则化牛顿下降方向的计算公式为：
33、 （6），
34、公式中，rg( λ)表示正则化牛顿下降方向；
35、s4.8.对待求物理变量进行迭代，并令，返回步骤s4.2，
36、迭代公式为：
37、 （7），
38、公式中， t( λ)表示正则化牛顿下降方向上的搜索步长；
39、s4.9.将最终变量 λ*代入s4.3的形函数公式中，计算得到局部最大熵插值函数，根据局部最大熵插值函数计算其导数，局部最大熵插值函数的导数的计算公式为：
40、 （8），
41、其中，（9），
42、（10），
43、公式中，r*表示原目标函数的梯度矢量，j*表示原目标函数的黑塞矩阵。
44、优选地，所述s4.8中的迭代过程，正则化牛顿下降方向上的搜索步长 t( λ)初步设置为1，所述s4.8对待求物理变量进行迭代后，计算迭代前后的目标函数值，若迭代后的目标函数值小于迭代前的目标函数值，则直接返回步骤s4.2，若迭代后的目标函数值不小于迭代前的目标函数值，通过精确线搜索算法计算正则化牛顿下降方向上的最优步长 t( λ)，并按照迭代公式重新进行迭代，重新迭代后返回步骤s4.2。
45、优选地，所述s4.1中待求物理变量的初始值 λ0通过以下方式获得：
46、s4.1.1.判断问题域所要解决的问题是三维问题还是二维问题，若是三维问题，在节点集中选择距离物质点x最近的4个节点组成一个四面体，计算物质点x关于四面体的四个顶点的重心坐标，若是二维问题，在节点集中选择距离物质点x最近的3个点组成一个三角形，计算物质点x关于三角形的三个顶点的重心坐标；
47、s4.1.2.令物质点的形函数等于相应的重心坐标的分量，分别代入形函数的计算公式中，整理后得到一组关于变量 λ的线性方程组，求解后得到变量的初始值 λ0。
48、优选地，所述s4.1.2求解得到变量的初始值 λ0后，对各节点的配分函数值进行求和得到 z( λ0)，判断该值是否小于最小阈值或大于最大阈值，若是，则采用nelder-mead算法对s4.1.1中获得的重心坐标进行修正，修正方式如下：
49、s4.1.2.1.在 λ空间内选取一个与 λ0关联的简单形，设置容差 ε和最大比较次数nmax，设置初始操作次数nfunk设为0；
50、s4.1.2.2.根据目标函数值对简单形顶点进行排序，将最差的顶点标记为w，次差的顶点标记为s，最好的顶点标记为b，计算最差顶点w处变量的目标函数值 f(w)与最好顶点b处变量的目标函数值 f(b)的差，若 f(w)- f(b)＞ ε或nfunk>nmax，则输出修正后的重心坐标，否则进入s4.1.2.3开始修正；
51、s4.1.2.3.抛弃最差顶点w，找到剩余顶点所组成的面或线的中心m，找到最差顶点w关于中心m的对称点r，操作次数nfunk加2；
52、s4.1.2.4.比较对称点处变量的目标函数值 f(r)与最差顶点w处变量的目标函数值 f(w)，若 f(r)＜ f(w)，则进入s4.1.2.5，否则进入s4.1.2.7；
53、s4.1.2.5.比较对称点处变量的目标函数值 f(r)与最好顶点b处变量的目标函数值 f(b)的大小，若 f(r)≤ f(b)，则在重心m与对称点r的延长线上确定e点，e点的计算方式为e=m+2（r-m），进一步比较e点处变量的目标函数值 f(e)与对称点处变量的目标函数 f(r)的大小，若 f(e)＜ f(r)，用e点替换最差的顶点w，返回s4.1.2.2，若 f(e)≥ f(r)，用对称点r代替最差的顶点w，返回s4.1.2.2；若 f(r)＞ f(b)，则进入s4.1.2.6；
54、s4.1.2.6.比较对称点处变量的目标函数 f(b)与次差顶点s处变量的目标函数 f(s)，若 f(r)≤ f(s)，用对称点r代替最差顶点w，操作次数nfunk减1，返回s4.1.2.2，若 f(r)＞ f(s)，进入s4.1.2.7；
55、s4.1.2.7.用对称点r代替最差顶点w，取重心m与对称点r的连线的中点c；
56、s4.1.2.8.进一步比较c点处变量的目标函数值 f(c)与最差顶点w处变量的目标函数值 f(w)，若 f(c)＞ f(w)，用c点代替最差顶点w，返回s4.1.2.2，若 f(c)≤ f(w)，令简单形向最好顶点b收缩，输出收缩后的顶点，操作次数nfunk加3，返回s4.1.2.2。
57、优选地，所述s4.2中的帕雷托最优化参数通过以下公式计算：
58、 （11），
59、公式中， h为节点集构成的凸集的局部特征长度，取最小节点间距， γ为无量纲参数；
60、所述s8中邻域节点集更新后，均重新计算帕雷托最优化参数。
61、优选地，所述s5的具体步骤包括：
62、s5.1.更新节点集中质量，公式为：
63、（12），
64、公式中，， n a表示局部最大熵插值函数， m表示节点集中质量， a表示节点索引号， n表示 n时刻， p表示物质点索引号， n h(x a, n)表示节点x a的邻域物质点集；
65、s5.2.更新下一个节点的节点坐标，表示为：
66、（13），
67、公式中，l a, n表示索引号为 a的节点在 n时刻的线性动量，f a, n表示索引号为 a的节点在 n时刻的节点力；
68、s5.3.基于下一个节点的节点坐标以及最大熵插值函数的导数计算局部变形梯度增量，计算公式为：
69、（14），
70、公式中，f p, n→n+1表示从 n时刻到下一时刻的局部变形梯度增量， n h(x p, n)表示每个物质点x p的邻域节点集；
71、s5.4.计算下一时刻物质点的局部变形梯度，计算公式为：
72、（15），
73、公式中，f p, n表示 n时刻的局部变形梯度，f p, n+1表示下一时刻的局部变形梯度。
74、优选地，所述s7中的更新的节点力包括节点内力和节点外力，节点力的计算公式为：
75、（16），
76、公式中，f int a,n+1表示索引号为 a的节点下一时刻的节点内力，f ext a,n+1表示索引号为 a的节点下一时刻的节点外力，f a,n+1表示索引号为 a的节点下一时刻的节点力；
77、所述s7中更新线性动量的计算公式为：
78、（17），
79、公式中，l a, n表示索引号为 a的节点在 n时刻的线性动量，l a, n+1表示索引号为 a的节点在下一时刻的线性动量。
80、优选地，所述s8中更新物质点坐标的公式为：
81、（18），
82、公式中，x a, n+1表示下一时刻的节点坐标，x p, n+1表示下一时刻的物质点坐标， n h(x p, n)表示物质点x p的邻域节点集；
83、更新下一个时刻邻域节点集 n h(x p, n+1)，具体通过以下方式获得：给定搜索半径，确定每个物质点在搜索半径内的节点为该物质点的邻域节点，将邻域节点输入为邻域节点集。遍历物质点的邻域节点集 n h(x p, n+1)，构建节点的邻域物质点集 n h(x a, n+1)。
84、本发明还涉及一种基于无网格法插值技术的冲击大变形过程模拟方法的应用，其应用于模拟子弹撞击靶板的冲击大变形过程。
85、采用本发明提供的技术方案，与现有技术相比，具有如下有益效果：
86、1.本发明涉及的基于无网格法插值技术的冲击大变形过程模拟方法采用正则化牛顿法对问题域中的待求物理变量 λ进行迭代更新，修复了原目标函数黑塞矩阵行列式值为零的区域，修复后在空间内任意一点皆可构造正则化牛顿方向进行求解，不会造成牛顿法失去作用的问题。
87、2.本发明涉及的基于无网格法插值技术的冲击大变形过程模拟方法在采用正则化牛顿法对问题域中的待求物理变量 λ进行迭代更新时，对于三维问题，在节点集中选择距离物质点x最近的4个节点组成一个四面体，计算物质点x关于四面体的四个顶点的重心坐标，对于二维问题，在节点集中选择距离物质点x最近的3个点组成一个三角形，计算物质点x关于三角形的三个顶点的重心坐标，再令物质点的形函数等于相应的中心坐标的分量，分别代入形函数的计算公式中，整理后得到一组关于变量 λ的线性方程组，求解后得到变量的初始值 λ0，该方法求出的初始值位于最优解的附近，经数次迭代就能迅速收敛，增加迭代效率。
88、3. 本发明涉及的基于无网格法插值技术的冲击大变形过程模拟方法求得初始值后，对于某些极端情况，使用nelder-mead方法把迭代点向极值点收缩，保证了迭代值落入极端情况时计算不会崩溃，即保证了整个算法的稳定性。
1.一种基于无网格法插值技术的冲击大变形过程模拟方法，其特征在于：包括以下步骤：
2.根据权利要求1所述的基于无网格法插值技术的冲击大变形过程模拟方法，其特征在于：所述s4中采用正则化牛顿法对问题域中的变量λ进行迭代求解的具体步骤为：
3.根据权利要求2所述的基于无网格法插值技术的冲击大变形过程模拟方法，其特征在于：所述s4.8中的迭代过程，正则化牛顿下降方向上的搜索步长t(λ)初步设置为1，所述s4.8对待求物理变量进行迭代后，计算迭代前后的目标函数值，若迭代后的目标函数值小于迭代前的目标函数值，则直接返回步骤s4.2，若迭代后的目标函数值不小于迭代前的目标函数值，通过精确线搜索算法计算正则化牛顿下降方向上的最优步长t(λ)，并按照迭代公式重新进行迭代，重新迭代后返回步骤s4.2。
4.根据权利要求2所述的基于无网格法插值技术的冲击大变形过程模拟方法，其特征在于：所述s4.1中待求物理变量的初始值λ0通过以下方式获得：
5.根据权利要求4所述的基于无网格法插值技术的冲击大变形过程模拟方法，其特征在于：所述s4.1.2求解得到变量的初始值λ0后，对各节点的配分函数值进行求和得到z(λ0)，判断该值是否小于最小阈值或大于最大阈值，若是，则采用nelder-mead算法对s4.1.1中获得的重心坐标进行修正，修正方式如下：
6.根据权利要求2所述的基于无网格法插值技术的冲击大变形过程模拟方法，其特征在于：所述s4.2中的帕雷托最优化参数通过以下公式计算：
7.根据权利要求2所述的基于无网格法插值技术的冲击大变形过程模拟方法，其特征在于：所述s5的具体步骤包括：
8.根据权利要求7所述的基于无网格法插值技术的冲击大变形过程模拟方法，其特征在于：所述s7中的更新的节点力包括节点内力和节点外力，节点力的计算公式为：
9.根据权利要求7所述的基于无网格法插值技术的冲击大变形过程模拟方法，其特征在于：所述s8中更新物质点坐标的公式为：
10.一种根据权利要求 1 所述的基于无网格法插值技术的冲击大变形过程模拟方法的应用，其应用于模拟子弹撞击靶板的冲击大变形过程。
