CN106446432A - 一种求解材料大变形的最优输运无网格方法 - Google Patents

Abstract

本发明涉及一种求解材料大变形的最优输运无网格方法(Optimal Transportation Meshfree，OTM)，为了高效稳定的求解极大变形、高速冲击及几何畸变、金属材料成型及多相耦合等问题。OTM方法采用物质点与节点对初始问题域进行离散，采用局部最大熵插值函数来构造连续的运动函数，避免了有限元方法处理极大变形网格畸变，无网格法无法直接添加Dirichlet边界条件，以及计算不收敛等问题。同时由于插值与积分在不同的离散点进行，提供了有效的无网格数值积分模式，解决了拉应力不稳定性的问题。OTM方法作为增量更新拉格朗日方法，质量守恒自动满足，无需求解。通过采用最优运输理论进行时间离散，保证了离散系统的动量守恒以及辛守恒，极大的提高了运算的效率与精度。

Description

一种求解材料大变形的最优输运无网格方法

技术领域

本发明涉及一种连续介质力学问题中求解材料动态响应的方法，特别是涉及一种连续介质力学问题中求解材料动态响应的无网格方法，属于计算力学与数值仿真技术领域，主要应用于求解极度大变形、动态裂纹扩展、高速冲击及几何畸变、材料裂变、金属材料成型及多相变等问题。

背景技术

在采用有限元方法处理现代工程系统中包括：极度大变形问题；动态裂纹扩展问题；高速冲击及几何畸变问题；材料裂变问题；金属材料成型问题及多相变问题等在内的复杂的物理现象时，由于巨大的网格畸变或单元分裂造成有限元求解困难甚至导致求解失败，为了解决这些问题，往往在有限元计算中不断地进行网格重新划分，然而，这样不但大大地增加了计算时间，而且对于有些问题仅仅重新划分网格并不能完全得以解决。有限元的这些缺点源自于它使用了事先定义好的网格，要想彻底解决有限元所面临的这些问题，就应避免应用固定的网格。因此近年来无网格法的思想被提了出来，并得到了迅速的发展，无网格法是在对于一个问题域建立离散的系统方程时不用事先定义好的网格的一种数值方法，与有限元法相比，无网格法的优点包括：不需要网格；容易构造高阶形状函数；容易进行自适应分析。无网格方法的这些优点使其具有相当的潜力解决上述提及的传统有限元法所不能处理的极度大变形、动态裂纹扩展、高速冲击及几何畸变、材料裂变、金属材料成型及多相变等问题。目前，根据无网格法所使用的计算模型不同，无网格法可分为3大类：a)基于配点(collocation)的无网格法；b)基于积分弱式(weak form)的无网格法；c)基于积分弱式-配点结合的无网格法。

基于配点的无网格法是指直接离散微分方程(或偏微分方程)的一类数值方法，这类方法包括：无网格vortex法(参见Bernard P S.A deterministic vortex sheet methodfor boundary lalyer flow.Journal of Computation of Physics,1995,117:132-145)，基于任意网格的差分法(参见Liszka T.The Finite Difference Method for ArbitraryIrregular Meshes,a Variational Approach to Applied Mechanics Problems.Paris:GAMNI 2,1980.227-235)，有限点法(参见Onate E.A finite point method incomputational mechanics.applications to convective transport and fluidflow.Int J Numer Methods Engrg.1996.39:3839-3867)，hp-云法(参见Armando D C.Hpclouds-a，meshless method to solve boundary value problems.TIC AM Report 95-05.University of Texas at Austin.1995)和各种形式的无网格配点法。这类方法具有形式简单，计算效率高，不需要积分网格的优点，然而这类方法的缺点也十分明显：这类方法精度较低，稳定性差，特别对于具有Neumann(导数)边界条件的偏微分方程的求解就更加困难，而这种边界条件在固体力学问题中大量存在(如应力边界条件就是Neumann条件)，造成这一问题的主要原因是基于配点的无网格法在处理应力边界条件时不能很好地控制误差，并容易造成求解的不稳定性，因此，到目前为止，基于配点的无网格法仍主要用于流体问题的求解中，极大的限制了其应用范围。

基于积分弱式的无网格方法，此类无网格方法中获得离散方程一般应用Galerkin方法，通过对原控制方程的弱形式实施Galerkin过程，然后应用无网格形状函数进行离散。这类方法包括：自由单元Galerkin法(参见Belytschko T.Element-free Galerkinmethods.Int J Numer Methods Engrg,1994,37:229-256)，无网格局部伽辽金法(参见Atluri S N.A new meshless local Petrov-Galerkin approach in computationalmechanics.Computational Mechanics 1998.22,117-127)，无网格点插值法(参见Gu YT.Development of meshfree techniques for computational mechanics:Ph Dthesis.National University of Singapore,2002)，边界点法(参见Chati M K.Theboundary node method for three-dimensional problem s in potential theory.IntJ Numer Methods Engrg,2000,47:1523-1547)。这类方法具有稳定性高和精度好的优点，并已成功地应用在固体力学问题的分析中，然而基于弱式的无网格方法缺点也十分明显，它们普遍存在数值积分困难和本质边界条件难处理的问题，例如：全域弱式法中需要用全域背景网格进行数值积分不能完全摆脱对网格的依赖、形函数不满足Kronecker-delta属性导致应用本质边界条件非常困难(如自由单元Galerkin法)；局部弱式法中因权函数和试探函数取自不同的空间和不对称的边界积分引起“刚度阵”带状但不对称，这将增加求解的难度和计算量，尤其是对于一些具有复杂边界的问题，局部积分变得更难处理(如无网格局部伽辽金法)；计算效率低形状函数难以满足全域相容性条件将影响收敛性和精度(如无网格点插值法)。

基于积分弱式-配点结合的无网格方法，是新近提出的结合了前述基于配点的无网格法和基于积分弱式的无网格方法的特性无网格方法(参见Gu Y T.A meshfree weak-strong form mthod for time dependent problems.Computational Mechanics,2005,35,134-145和Liu G R.A meshfree method:meshfree weak-strong form method,for 2-D solids.Computational Mechanics,2003,33:2-14)。这种方法应用了配点和局部积分弱形式相结合的思想，使得其具有兼备配点法形式简单和计算效率高，及无网格弱式法精度高稳定性好的潜能，不过由于才提出不久，还有许多技术问题有待进一步解决，例如如何进一步提高稳定性，如何提高h收敛性，应用其解决复杂的问题等，还需要更多的研究工作。

综上所述，目前主流的无网格方法存在的最主要缺点是：

第一，不满足克罗内克属性。目前绝大部分无网格方法的插值函数不满足克罗内克属性，这将使得计算精度较低，同时由于离散域边界上的节点形函数值不为零，使得施加本质边界条件非常困难，需要经过复杂的特殊处理。

第二，“弱”形式下的等效数值积分困难。在缺少网格的条件下，在对无网格伽辽金法进行“弱”形式下的等效积分时会有较大困难，大量研究文献表明，目前无网格方法中主要采用的几种积分方案，如：节点积分方案、应力点积分方案、背景网格积分方案，它们都不同程度的面临着处理拉应力不稳定性、数值积分误差等困难。

第三，计算量大，效率低。无网格方法使用了高阶插值函数，提高了插值的连续性，但同时也极大地提高了计算量，求解域内的每个点的插值系数都需要对相应节点耦合矩阵求逆，光滑性越高耦合矩阵的阶数就越高，相应的求逆就越耗费计算量。

如前所述，无网格方法不需要网格(或对网格依赖性很小)、容易构造高阶形状函数、容易进行自适应分析的优点使其具有相当的潜力解决现代工程系统中极度大变形、动态裂纹扩展、高速冲击及几何畸变、材料裂变、金属材料成型及多相变等问题，然而目前已发展的无网格方法各有其优缺点，普遍存在边界条件难施加、数值积分困难及计算效率低的缺点，限制了其应用，仍需发展更加有效理想的无网格方法。

发明内容

本发明技术解决方案：克服现有技术的不足，提供一种求解材料大变形的最优输运无网格方法(Optimal Transportation Meshfree method，下文简称为OTM方法)，通过物质点和节点相结合的方式对问题域进行离散，其形状函数(局部最大熵插值函数)在边界上满足克罗内克属性，克服了绝大部分无网格法在处理位移边界条件时的困难，局部最大熵插值函数值及其导数可以通过牛顿迭代法高效求解，求解效率与邻域中节点的个数无关，避免了求解一般无网格插值函数的高计算量，提高了计算效率。同时其最优输运时间积分方案保证了时间离散下系统的质量与动量守恒以及辛守恒，由于插值与积分在不同的离散点进行，解决了拉应力不稳定性的问题。

本发明的技术解决方案：一种求解材料大变形的最优输运无网格方法，物质点与节点相结合的几何离散方案，最优输运理论时间离散方案，形状函数为局部最大熵插值函数，其在边界满足克罗内克属性，实现步骤如下：

第一步，设Ω表示d维的连续介质问题域，即几何模型，将几何模型Ω离散为一组物质点集{x_p,k,p＝1,2,…,m；k＝0,1,…,N}和一组节点集{x_a,k,a＝1,2,…,n；k＝0,1,…,N}，其中下标p和a代表物质点与节点，m和n为物质点与节点的个数，k代表时间步；

第二步，在初始时刻k＝0时，利用第一步离散得到的节点集和物质点集，初始化节点坐标x_a,k、节点形函数N_a(x_p,k)、节点形函数导数节点力f_k、节点线性动量矩阵l_k、节点质量矩阵M_k，初始化物质点坐标x_p,k、物质点邻域、物质点体积v_p,k、物质点密度ρ_p,k、物质点变形梯度F_p,k；

第三步，在第二步中初始化了k时刻节点集数据基础上，采用最优输运理论进行时间离散，对节点进行k→k+1时刻的瞬态分析，显式计算节点在k+1时刻的坐标

第四步，根据第三步得到的节点坐标x_k+1，利用局部最大熵插值函数得到物质点在k+1时刻的局部变形数据，包括：物质点增量变形梯度和物质点变形梯度F_p,k+1＝F_p,k→k+1οF_p,k；

第五步，根据第四步得到的物质点在k+1时刻的物质点增量变形梯度和物质点变形梯度，更新k+1时刻的节点力f_k+1，节点动量矩阵l_k+1、节点质量矩阵M_k+1；

第六步，根据第五步得到物质点从k→k+1时刻运动和变形数据，更新物质点的坐标x_p,k+1＝x_a,k+1N_a(x_p,k)、体积v_p,k+1＝det(F_p,k→k+1)v_p,k与密度其中m_p为物质点质量；

第七步，根据第六步计算得到的物质点坐标x_p,k+1、体积v_p,k+1与密度ρ_p,k+1，重新计算物质点邻域，并更新邻域内节点的插值函数值N_a(x_p,k+1)及导数值

第八步，在完成第七步后，即代表完成了物质点和节点在一个时间步内的动态分析，判断当前的时间步k，如果k＝N，代表已计算至最后一个时间步，此时退出计算，如果k≠N则转入第三步。

本发明与现有技术相比的优点在于：

(1)本发明采用物质点与节点相结合的方式对几何模型进行离散，计算过程不涉及网格，与传统的有限元方法相比完全避免了处理极度大变形、动态裂纹扩展、高速冲击及几何畸变、材料裂变、金属材料成型及多相变等问题时，由于固定网格所带来的问题，同时由于插值在节点进行而积分在物质点进行，解决了绝大部分无网格方法中存在的拉应力不稳定性；

(2)本发明采用局部最大熵无网格插值函数，在边界上满足克罗内克属性，在计算中可以像有限元方法一样直接添加位移边界条件，解决了绝大部分无网格法在处理位移边界条件时需要进行复杂计算的缺陷，同时，局部最大熵插值函数值及其导数通过牛顿迭代法高效求解，求解效率与邻域中节点的个数无关，避免了求解一般无网格插值函数的高计算量，提高了计算效率；

(3)本发明采用最优输运理论进行时间离散，保证了时间离散下系统的动量守恒与辛守恒自动满足无需额外计算，克服了欧拉方法中需要通过求解复杂泊松方程来达到质量守恒的高计算成本问题。

附图说明

图1为本发明一种求解材料大变形的最优输运无网格方法流程图；

图2为本发明应用于冲击速度从100m/s到10km/s的金属靶板冲击试验仿真，其中(a)代表弹道冲击速度为100m/s到600m/s的试验结果(左)和仿真结果(右)，(b)代表弹道冲击速度为1km/s到3km/s的高速冲击仿真结果，(c)代表弹道冲击速度为5km/s到10km/s的超高速冲击仿真结果；

图3为本发明通过物质点与节点相结合的方式对几何模型进行空间离散的示意图；

图4为本发明所采用的拉格朗日动力学时间离散方案；

图5为本发明离散后的几何模型从k时刻变形到k+1时刻的示意图；

图6为本发明所采用的局部最大熵插值函数在取不同γ值时，形状函数及其导数性质的二维例子示意图，其中，第一行图形为形状函数，第二行为形状函数在y方向导数，第三行为形状函数在x方向上的导数。

具体实施方式

下面结合附图，对本发明一种求解材料大变形的最优输运无网格方法的技术方案做进一步说明。

如图1所示，本发明具体实现如下：

(1)几何模型离散，本发明适用于各种材料模型，具体计算过程中用到哪种材料需根据用户实际情况而定，在进行材料动态响应分析时(如图2所示的金属靶板冲击试验)，设Ω表示d维的连续介质问题域(即几何模型)将几何模型Ω离散为一组物质点集{x_p,k,p＝1,2,…,m；k＝0,1,…,N}和一组节点集{x_a,k,a＝1,2,…,n；k＝0,1,…,N}。如图3所示，空心点代表节点：x_a,k(a为离散域中节点索引号，代表第几个节点，k代表第几个时间步，时间步长及步数由用户控制)，实心点代表物质点：x_p,k(p为离散域中物质点索引号，代表第几个物质点，k代表第几个时间步)，物质点与节点的初始位置可由用户根据不同的算法来确定，比如随机插入物质点与节点，以距离物质点最近的d+1个节点作为其邻域(d为分析问题的维度)，如图3中圆圈包含的3个节点即为箭头所指物质点的邻域。在本发明的具体实现中，为了保证计算的精度，采用三角形单元(二维情况)或者四边形单元(三维情况)对几何模型进行网格划分生成初始节点，每个单元的形心取为物质点。该网格只用于初始化节点与物质点的位置，将不会出现在计算过程中。

(2)初始化计算参数，在初始时刻k＝0时，利用第一步离散得到的节点集和物质点集，初始化节点坐标x_a,k、节点形函数N_a(x_p,k ⁾、节点形函数导数节点力f_k、节点线性动量矩阵l_k、节点质量矩阵M_k，初始化物质点坐标x_p,k、物质点邻域、物质点体积v_p,k、物质点密度ρ_p,k、物质点变形梯度F_p,k；

(3)计算k+1时刻节点坐标，在第(2)步中初始化了k时刻节点集数据，本发明采用最优输运理论进行时间离散，对节点进行k→k+1时刻的瞬态分析，显式预估节点在下一个时间步长时的坐标：

在公式(1)中节点从时间步t_k到t_k+1的坐标变化涉及到对时间的离散，OTM方法中采用最优输运理论对时间进行离散。

Monge-Kantotovich最优质量输运问题可描述为：质量相等的物质给定两种不同的质量分布g₀(x)和g₁(x)，找到一种最优的映射关系ψ，使得将物质从按g₀(x)质量分布输运至按g₁(x)质量分布所需的成本最小，如公式(2)所示：

同时两种分布要满足质量相等的条件，即：

更一般地，可以将分布函数g₀和g₁用非负值分布ρ₀和ρ₁进行替换，ρ₁＝ψ_#ρ₀，同时将公式(2)中的欧式距离|x-ψ(x)|，采用更为一般的成本函数c(x,y)进行表示，得到最优输运成本问题求解的一般形式：

因此质量输运的成本大小取决于两种质量分布之间的c(x,ψ(x))距离，这个距离被称为Kantorovich-Wasserstein距离。

对于空间中的一般连续介质的运动求解问题，在时间范围[a,b]，物质的运动需要满足耦合控制方程(5)和(6)：

其中，ρ是物质的密度，v是速度场，σ为Cauchy应力张量，b为体积力，公式(5)表示的是质量守恒方程，公式(6)表示的是线性动量守恒方程。假定物质的总质量M是有限的：

由于物质在运动过程中质量是守恒的，其质量是一个常数值，因此，物质的质量流动将在一个有限范围内进行，并不会扩散到无限大的空间中即：

其中B_R是半径为R的球体，n为速度场的法向，由于质量守恒，所以可以得到质量M的变化率为0，即：

在物质运动过程中，给定初始时刻和结束时刻物质密度：

ρ(x,a)＝ρ_a(x) (10)

ρ(x,b)＝ρ_b(x) (11)

因此物质从a时刻到b时刻的运动问题，等同于将物质从a时刻按ρ_a(x)质量分布输运至b时刻按ρ_b(x)质量分布，所以，这是前面阐述的质量输运问题，可以采用最优输运理论进行描述与求解。在此基础上，Benamou和Brenier证明了公式(5)至(11)所描述的质量输运问题具有变分特性，因此公式(5)至(11)所描述的问题等效于：

表示对于一般连续介质运动，最优的输运方案是使得输运过程所需要的动能和体积能最小，约束条件为质量守恒方程(5)，其中：

因而对于一般连续介质的运动描述，可以表示为质量的最优输运成本：

式中T₂(ρ_a,ρ_b)表示Kantorovich-Wasserstein距离，即物质从ρ_a(x)输运至ρ_b(x)的最优输运距离，用符号表示：

在本发明中，最优输运定理进行时间离散的实现方式如下：

(3.1)采用最优输运理论得到半离散化时间积分形式

在实际中物质从ρ_a(x)输运至ρ_b(x)是个连续的过程，为了进行求解需要将[a,b]这个时间范围进行离散，通过求解每个时间步的输运过程，最后近似得到实际连续的输运过程，根据物质从ρ_a(x)输运至ρ_b(x)可以有许多不同的路径如图4所示，因此有不同的离散方式，公式(16)中最优输运理论给出是最优的输运路径，对公式(14)按公式(16)给出的最优路径进行时间离散，得到半离散形式的时间积分：

(3.2)计算每个时间步的Kantorovich-Wasserstein距离

公式(17)中第一项为离散化的动能，通过求解每一个时间步的质量守恒和动量守恒方程获得其值，最后将所有时间步计算得到的联系起来得到最优输运路径。

(3.3)计算每个时间步的体积能

公式(17)中第二项代表离散化的体积能，体积能的离散方式采用Center Difference的方式进行离散。

(3.4)显式计算节点在k+1时刻的坐标x_k+1

公式(17)直接是密度ρ的函数，通过4.2)和4.3)计算得到A_d在{ρ₁,…s,ρ_N-1}∈[L¹(Rⁿ；[0,∞))]^N-1的精确解，在[t_k,t_k+1]内该离散运动由增量传输映射实现从ρ_k至ρ_k+1的最优传输。

图5所示即为OTM方法中采用最优输运理论对计算邻域的拉格朗日函数进行时间离散，给定初始条件便可通过增量迭代的方式得到系统动力学参数的时间历程，其中包括节点的坐标，速度、加速度随时间的变化等信息，据此完成在k+1时刻的节点坐标x_k+1的更新。

(4)根据第(3)步计算得到的节点坐标x_k+1，此时节点就已产生一个位移，同时物质点产生变形，物质点的运动将通过节点的动力学信息动态插值获得，从而得到每个时间步长材料的局部变形，具体为：利用局部最大熵插值函数计算物质点上的局部增量变形梯度及物质点变形梯度F_p,k+1＝F_p,k→k+1°F_p,k来计算物质点的变形。

在构造局部最大熵插值函数时，设点集X为由n个节点组成的凸集，凸集X中每个节点x_a的形函数定义为N_a，该形函数在给定的抽样点x的值可通过最小化公式(19)求得：

同时对于需满足条件：

N_a(x)≥0,a∈[1,n] (20)

ΣN_a(x)＝1 (21)

ΣN_a(x)x_a＝x (22)

公式(20)至(22)保证了形函数的非负性以及零阶和一阶连续性，局部最大熵形函数是精确的线性插值，公式(19)的最小值唯一，其解即为局部最大熵形函数，如公式(23)所示：

其中，

在公式(23)中为多目标优化方法中的帕雷托最优参数(Pareto)，其大小决定了该多目标优化问题的最优解；通过求解最小化公式(19)的Dual形式得到λ^*(x)，进而获得Z(x,λ^*(x))的最小值，这种无约束凸优化问题可以通过Newton-Raphson迭代法和Nelder-Mead单纯形算法相结合高效稳定求解。局部最大熵形函数的精度会随着节点密度的增加而增加，其收敛性具有严格的数学证明。

在本发明中，局部最大熵插值函数的构造方法如下：

(4.1)定义维度无关参数γ

当计算所有节点的形函数时，采用与维度无关的参数γ来度量公式(23)中熵值项的局部性，图6所示即为通过调整γ的不同取值进而控制形函数的局部性。

(4.2)计算帕雷托最优化参数β

根据物质点大小h及γ大小，计算多目标优化帕雷托最优化参数β，公式为：

其中h为凸集X中的局部特征长度。

(4.3)计算参数λ^*

通过结合Newton-Raphson迭代法和Nelder-Mead单纯形算法求解最小化公式(19)的Dual形式得到λ^*(x)，进而获得Z(x,λ^*(x))的最小值。

(4.4)计算局部最大熵插值函数N_a(x)和导数

其中，

在OTM方法中，物质点运动通过节点的动力学信息插值获得，公式(27)中x即为物质点，x_a为节点，N_a(x)即代表节点x_a的形函数在物质点x的取值。

(5)根据第(4)步计算得到的物质点局部变形，通过材料的本构关系得到物质点的应力，并计算该增量运动引起的节点力_fk+1，同时更新节点的动量_lk+1与质量_Mk+1；

由更新后的变形梯度F_p,k+1根据材料本构关系得到局部应力P_iJ(F_p,k+1)，然后根据公式(30)计算节点力：

其中，x_p,k为t_k时刻物质了的坐标，为物质点t_k时刻的总变形梯度，为形函数导数，w_p,k为物质点体积。更新节点力后，根据公式(31)更新节点动量矩阵和公式(32)更新节点质量矩阵：

(6)根据第(5)步得到物质点从k→k+1时刻运动和变形数据，更新物质点的坐标x_p,k+1＝x_a,k+1N_a(x_p,k)、体积v_p,k+1＝det(F_p,k→k+1)v_p,k与密度

(7)根据第(6)步计算得到的物质点坐标x_p,k+1、体积v_p,k+1与密度ρ_p,k+1，采用cellarray算法在物质点邻域搜索范围内进行搜索判断新增加的节点，并将符合作为物质点邻域要求的节点增加至原邻域中，实现对每个物质点每个时间步邻域范围的动态增量更新，并计算邻域内节点的插值函数值N_a(x_p,k+1)及导数值

(8)在完成第(7)步后，即代表完成了物质点和节点在一个时间步内的动态分析，判断当前的时间步k，如果k＝N，代表已计算至最后一个时间步，此时退出计算，如果k≠N则转入第(3)步。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (4)

1.一种求解材料大变形的最优输运无网格方法，其特征在于包括以下步骤：

第一步，设Ω表示d维的连续介质问题域，即几何模型，将几何模型Ω离散为一组物质点集{x_p,k,p＝1,2,…,m；k＝0,1,…,N}和一组节点集{x_a,k,a＝1,2,…,n；k＝0,1,…,N}，其中下标p和a代表物质点与节点，m和n为物质点与节点的个数，k代表时间步；

第二步，在初始时刻k＝0时，利用第一步离散得到的节点集和物质点集，初始化节点坐标x_a,k、节点形函数N_a(x_p,k)、节点形函数导数节点力f_k、节点线性动量矩阵l_k、节点质量矩阵M_k，初始化物质点坐标x_p,k、物质点邻域、物质点体积v_p,k、物质点密度ρ_p,k、物质点变形梯度F_p,k；

第三步，在第二步中初始化了k时刻节点集数据基础上，采用最优输运理论进行时间离散，对节点进行k→k+1时刻的瞬态分析，显式计算节点在k+1时刻的坐标

第四步，根据第三步得到的节点坐标x_k+1，利用局部最大熵插值函数得到物质点在k+1时刻的局部变形数据，包括：物质点增量变形梯度和物质点变形梯度F_p,k+1＝F_p,k→k+1оF_p,k；

第五步，根据第四步得到的物质点在k+1时刻的物质点增量变形梯度和物质点变形梯度，更新k+1时刻的节点力f_k+1，节点动量矩阵l_k+1、节点质量矩阵M_k+1；

第六步，根据第五步得到物质点从k→k+1时刻运动和变形数据，更新物质点的坐标x_p,k+1＝x_a,k+1N_a(x_p,k)、体积v_p,k+1＝det(F_p,k→k+1)v_p,k与密度其中m_p为物质点质量；

第七步，根据第六步计算得到的物质点坐标x_p,k+1、体积v_p,k+1与密度ρ_p,k+1，重新计算物质点邻域，并更新邻域内节点的插值函数值N_a(x_p,k+1)及导数值

第八步，在完成第七步后，即代表完成了物质点和节点在一个时间步内的动态分析，判断当前的时间步k，如果k＝N，代表已计算至最后一个时间步，此时退出计算，如果k≠N则转入第三步。

2.根据权利要求1所述的一种求解材料大变形的最优输运无网格方法，其特征在于：所述第一步中，采用节点与物质点对问题域进行离散具体如下：

步骤1：二维问题采用三角形单元/三维问题采用四边形单元对几何模型进行网格划分生成初始节点，得到离散的节点；

步骤2)：根据步骤1)划分的单元，取每个单元的形心取为物质点，得到离散的物质点；

步骤3)：由步骤1)和步骤2)所得到的节点和物质点，将材料的物理信息存储在物质点上，包括变形、应力、材料内部参数，将材料的动力学信息存储在节点上，包括位移、速度、加速度与温度，计算过程不涉及网格，插值在节点进行，积分在物质点进行。

3.根据权利要求1所述的一种求解材料大变形的最优输运无网格方法，其特征在于：所述第三步中最优输运理论进行时间离散具体如下：

步骤1)：采用最优输运理论得到半离散化时间积分形式；

步骤2)：通过求解质量守恒和动量守恒方程，得到半离散化时间积分形式中每个时间步的Kantorovich-Wasserstein距离；

步骤3)：计算每个时间步的体积能，体积能采用Center Difference离散得到；

步骤4)：由步骤2)和步骤3)精确求解步骤1)中的半离散化时间积分，显式计算节点在k+1时刻的坐标x_k+1。

4.根据权利要求1所述的一种求解材料大变形的最优输运无网格方法，其特征在于：所述第四步中采用的局部最大熵插值函数的构造方法如下：

步骤1)：由用户输入维度无关参数γ大小，取值范围为0.1到6.8；

步骤2)：根据物质点大小h及γ大小，计算多目标优化帕雷托最优化参数β，计算公式为

步骤3)：通过Newton-Raphson迭代法和Nelder-Mead单纯形算法求解参数λ^*，计算公式为其中参数λ代表拉格朗日乘子，参数λ^*代表通过最小化log Z(x_p,k,λ)获得的拉格朗日乘子；

步骤4)：计算局部最大熵插值函数 以及导数其中J(x_a,k,λ^*)是Hessian目标函数；

步骤5)：利用局部最大熵插值函数计算物质点上的局部增量变形梯度 及物质点变形梯度F_p,k+1＝F_p,k→k+1.F_p,k来计算物质点的运动和变形。