求解强热流固耦合问题的并行无网格方法及系统与流程

本发明涉及计算力学与数值仿真领域，特别是涉及一种求解强热流固耦合问题的并行无网格方法及系统。

背景技术：

在包括制造业、新能源产业等新兴产业在内的大部分工业应用中，材料通常都要经受强热流固耦合问题，尤其近来增材制造、三维打印等技术的迅速发展使得强热流固耦合问题成为了全球研究的重点之一。这些强热流固耦合问题通常包括了：极度大变形问题、大温度梯度问题、高加热或冷却速率问题、高速冲击及几何畸变问题、材料裂变问题、金属材料成型问题、多相变问题等复杂的物理现象。在采用传统的拉格朗日有限元方法处理这些问题时，由于巨大的网格畸变或单元分裂造成有限元求解困难甚至导致求解失败，为了解决这些问题，有限元计算中往往需要不断地进行网格重新划分，然而，这样不仅大大地增加了计算时间，而且对于有些问题仅仅重新划分网格并不能完全得以解决。另外一个解决该问题的传统方法为ale(arbitrarylangrangian-eulerianformulation)方法。虽然该方法兼具拉格朗日方法和欧拉方法的优点，它在解决大变形问题和非受限多相流问题时仍十分受限，且同样需要进行网格重新划分。

近年来基于拉格朗日方法的无网格法得到了迅速的发展。与有限元法相比，无网格法的优点包括：不需要网格；容易构造高阶形状函数；容易进行自适应分析。无网格方法的这些优点使其具有相当的潜力解决上述提及的传统有限元法所不能处理的大温度梯度、高加热或冷却速率、高速冲击及几何畸变、材料裂变、金属材料成型及多相变等问题。目前，无网格法可粗略分为3大类：a)基于配点(collocation)的无网格法；b)基于积分弱式(weakform)的无网格法；c)基于积分弱式-配点结合的无网格法。

基于配点的无网格法采用直接离散微分方程(或偏微分方程)的方法，这类方法包括：无网格vortex法、基于任意网格的差分法、有限点法、hp-云法和各种形式的无网格配点法。这类方法形式简单，计算效率高，不需要积分网格，然而精度较低，稳定性差，特别在处理neumann(导数)边界条件时误差较大，求解不稳定。因此，到目前为止，基于配点的无网格法仍主要用于流体问题的求解中，极大的限制了其应用范围。

基于积分弱式的无网格方法一般应用galerkin方法，通过对原控制方程的弱形式实施galerkin过程，然后应用无网格形状函数进行离散。这类方法包括：自由单元galerkin法、无网格局部伽辽金法、无网格点插值法以及边界点法。这类方法稳定性高、精度较好，却难以处理本质边界条件，尤其是对于一些具有复杂边界的问题，局部积分变得更难处理(如无网格局部伽辽金法)。

基于积分弱式-配点结合的无网格方法则结合了前述两种方法的特性。这种方法应用了配点和局部积分弱形式相结合的思想，兼具两种方法的优点。然而有许多技术问题有待进一步解决，例如如何进一步提高稳定性，如何提高h收敛性，应用其解决复杂的问题等，还需要更多的研究工作。

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

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

第二，“弱”形式下的等效数值积分困难，目前无网格方法中主要采用的几种积分方案，如：节点积分方案、应力点积分方案、背景网格积分方案，它们都不同程度的面临着处理拉应力不稳定性、数值积分误差等困难。

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

第四，缺乏强热流固耦合算法可同时预测大变形下的温度场与材料相变。

如前所述，无网格方法不需要网格(或对网格依赖性很小)、容易构造高阶形函数、容易进行自适应分析的优点使其具有相当的潜力解决现代工程系统中强热力学耦合问题中极度大变形、大温度梯度、高加热或冷却速率、高速冲击及几何畸变、材料裂变、金属材料成型及多相变等问题。然而，目前无网格法仍然主要用于求解力学问题，少有将无网格法应用于求解强热力学耦合问题的研究，也未有人提出基于无网格法实现求解包括能量耗散机制在内的非平衡热力学耦合问题的数学架构。因此，如何进一步发展无网格法来解决强热力学耦合问题，仍需要进行进一步研究。

技术实现要素：

本发明的目的是提供一种求解强热流固耦合问题的并行无网格方法及系统，提高求解强热力学耦合问题的运算的效率与精度。

为实现上述目的，本发明提供了如下方案：

一种求解强热流固耦合问题的并行无网格方法，包括：

获取d维的连续介质问题域；

将所述d维的连续介质问题域进行离散，得到物质点集和节点集；所述物质点集包括多个物质点；所述节点集包括多个节点；所述物质点用于表示某个点设定范围内的物质的物理信息；所述物理信息包括质量以及密度；所述节点用于表示某个点的动力学的信息；所述动力学的信息包括速度以及加速度；

获取所述d维的连续介质问题域的计算步数以及处理器的个数；

根据所述处理器的个数对所述物质点集进行划分，并将划分后的物质点集分别发送到各处理器中；

根据tk时刻每一所述处理器上各个物质点与所述节点集的坐标以及所述物质点的邻域范围确定每个所述物质点的邻域节点；

根据所述物质点的邻域节点确定所述处理器上所有物质点邻域所覆盖的空间范围；所述空间范围的上下边界值为一个长方体，所述空间范围为所述处理器计算域的rangebox；

根据所述处理器所有物质点邻域所覆盖的空间范围确定各个所述处理器相互重叠的范围；

根据各个所述处理器相互重叠的范围确定重叠节点；

给定tk时刻各个所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度，给定tk时刻的物质点坐标，以及采用局部最大熵插值函数计算tk时刻物质点的局部变形数据已知，所述局部形变数据包括物质点形函数取值、物质点形函数导数取值、物质点体积、物质点密度、物质点变形梯度和物质点温度；

根据所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度确定全局质量矩阵、全局节点力和全局加速度；

根据所述全局质量矩阵、全局节点力和全局加速度，采用最优输运理论进行时间离散，更新tk+1时刻的节点坐标和节点速度；

根据所述tk时刻的物质点的局部形变、所有节点在tk+1时刻的坐标，采用牛顿迭代法确定该处理器上tk+1时刻所述局部节点的热驱动力与局部温度刚度矩阵；

根据各处理器上所得的局部节点的热驱动力及局部温度刚度矩阵，组合成为全局节点上的热驱动力和全局温度刚度矩阵，并计算得到全局节点的温度场；

根据所述tk+1时刻所述节点的坐标确定更新后该处理器的rangebox，并根据所述处理器的rangebox确定各个所述处理器相互重叠的范围；

根据在tk+1时刻所有节点的坐标，更新物质点的坐标、体积和密度，并确定在tk+1时刻物质点的局部变形数据，所述局部变形数据包括物质点形函数取值、物质点形函数导数取值、物质点变形梯度和物质点温度；

根据在tk+1时刻所述物质点的局部变形数据确定tk+1时刻的节点的形变数据；所述节点的形变数据包括节点力、节点动量矩阵以及节点质量矩阵；

判断k+1是否等于所述计算步数；

若k+1等于所述计算步数，则根据所述物质点和节点从t0时刻到tk+1时刻形变数据和坐标变化以及温度场变化确定强热流固耦合问题的解，完成动态响应分析；

若k+1不等于所述计算步数，则令k＝k+1，并返回根据所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度确定全局质量矩阵、全局节点力和全局加速度的步骤。

一种求解强热流固耦合问题的并行无网格系统，包括：

第一获取模块，用于获取d维的连续介质问题域；

物质点集和节点集确定模块，用于将所述d维的连续介质问题域进行离散，得到物质点集和节点集；所述物质点集包括多个物质点；所述节点集包括多个节点；所述物质点用于表示某个点设定范围内的物质的物理信息；所述物理信息包括质量以及密度；所述节点用于表示某个点的动力学的信息；所述动力学的信息包括速度以及加速度；

第二获取模块，用于获取所述d维的连续介质问题域的计算步数以及处理器的个数；

划分模块，用于根据所述处理器的个数对所述物质点集进行划分，并将划分后的物质点集分别发送到各处理器中；

邻域节点确定模块，用于根据tk时刻每一所述处理器上各个物质点与所述节点集的坐标以及所述物质点的邻域范围确定每个所述物质点的邻域节点；

空间范围确定模块，用于根据所述物质点的邻域节点确定所述处理器上所有物质点邻域所覆盖的空间范围；所述空间范围的上下边界值为一个长方体，所述空间范围为所述处理器计算域的rangebox；

相互重叠的范围确定模块，用于根据所述处理器所有物质点邻域所覆盖的空间范围确定各个所述处理器相互重叠的范围；

重叠节点确定模块，用于根据各个所述处理器相互重叠的范围确定重叠节点；

初始化模块，用于初始化tk时刻的物质点坐标；设定k＝0；

tk时刻的物质点的局部形变数据确定模块，用于给定tk时刻各个所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度，给定tk时刻的物质点坐标，以及采用局部最大熵插值函数计算tk时刻物质点的局部变形数据，所述局部形变数据包括物质点形函数取值、物质点形函数导数取值、物质点体积、物质点密度、物质点变形梯度和物质点温度；

节点局部信息确定模块，用于根据所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度；

节点全局信息确定模块，用于根据所述全局质量矩阵、全局节点力和全局加速度，采用最优输运理论进行时间离散，更新tk+1时刻的节点坐标和节点速度；

tk+1时刻的节点坐标更新模块，用于根据所述tk时刻的物质点的局部形变、所有节点在tk+1时刻的坐标，采用牛顿迭代法确定所述处理器上tk+1时刻所述局部节点的热驱动力与局部温度刚度矩阵；

温度场确定模块，根据各处理器上所得的局部节点的热驱动力及局部温度刚度矩阵，组合成为全局节点上的热驱动力和全局温度刚度矩阵，并计算得到全局节点的温度场；

tk+1时刻rangebox更新模块，用于根据tk+1时刻所述节点的坐标确定更新后该处理器的rangebox，并根据所述处理器的rangebox确定各个所述处理器相互重叠的范围；

tk+1时刻物质点的局部变形数据确定模块，用于根据在tk+1时刻所有节点的坐标，更新物质点的坐标、体积、密度，并确定在tk+1时刻物质点的局部变形数据，所述局部形变数据包括物质点形函数取值、物质点形函数导数取值、物质点变形梯度和物质点温度；

tk+1时刻的节点的形变数据确定模块，用于根据在tk+1时刻所述物质点的局部变形数据确定tk+1时刻的节点的形变数据；所述节点的形变数据包括节点力、节点动量矩阵以及节点质量矩阵；

第一判断模块，用于判断k+1是否等于所述计算步数；

强热流固耦合问题求解完成模块，用于若k+1等于所述计算步数，则根据所述物质点和节点从t0时刻到tk+1时刻形变数据和坐标变化以及温度场变化确定强热流固耦合问题的解，完成动态响应分析；

空间范围更新模块，用于若k+1不等于所述计算步数，则令k＝k+1，并返回根据所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度确定全局质量矩阵、全局节点力和全局加速度的步骤。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明所提供的一种求解强热流固耦合问题的并行无网格方法及系统，通过采用物质点与节点对初始问题域进行离散，采用局部最大熵插值函数来构造连续的形函数，避免了有限元方法处理极大变形时网格畸变、无网格法无法直接添加dirichlet边界条件，以及计算不收敛等问题。本发明通过采用最优运输理论进行时间离散，保证了离散系统的动量守恒以及能量守恒，极大的提高了运算的效率与精度。本发明分布式并行化实现，更进一步提高了运算的效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明所提供的一种求解强热流固耦合问题的并行无网格方法流程示意图；

图2为本发明所提供的一种求解热力学耦合问题中材料大变形和温度场的无网格方法流程图；

图3为本发明所提供的一种求解热力学耦合问题中材料大变形和温度场的无网格方法分布式并行化实现过程流程图；

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

图5为本发明rangebox的定义方案示意图；

图6为本发明中求解能量平衡方程中对全局线性方程组刚度矩阵的分布式并行化示意图；

图7为本发明应用于求解极高速激光材料沉积制造技术的结果图；

图8为本发明所提供的拉格朗日动力学时间离散方案示意图；

图9为本发明离散后的连续介质问题域的几何模型从tk时刻变形到tk+1时刻的示意图；

图10为本发明所提供的局部最大熵插值函数在取不同γ值时，形状函数及其导数性质的二维例子示意图；

图11为本发明所提供的一种求解强热流固耦合问题的并行无网格系统流程示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种求解强热流固耦合问题的并行无网格方法及系统，提高求解强热力学耦合问题的运算的效率与精度。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明所提供的一种求解强热流固耦合问题的并行无网格方法流程示意图，如图1所示，本发明所提供的一种求解强热流固耦合问题的并行无网格方法，包括：

s101，获取d维的连续介质问题域；

s102，将所述d维的连续介质问题域进行离散，得到物质点集和节点集；所述物质点集包括多个物质点；所述节点集包括多个节点；所述物质点用于表示某个点设定范围内的物质的物理信息；所述物理信息包括质量以及密度；所述节点用于表示某个点的动力学的信息；所述动力学的信息包括速度以及加速度；s102具体包括：

根据所述d维的连续介质问题域，采用不同形状的单元进行网格划分，生成初始节点，得到离散的节点；不同形状的单元包括三角形单元和四面体单元；

根据所述单元的几何中心确定所述物质点；

将所述动力学的信息存储在所述节点上；

将所述物理信息存储在所述物质点上。

s103，获取所述d维的连续介质问题域的计算步数以及处理器的个数；

s104，根据所述处理器的个数对所述物质点集进行划分，并将划分后的物质点集分别发送到各处理器中；

s105，假定k＝0，根据tk时刻每一所述处理器上各个物质点与所述节点集的坐标以及所述物质点的邻域范围确定每个所述物质点的邻域节点；

s106，根据所述物质点的邻域节点确定所述处理器上所有物质点邻域所覆盖的空间范围；所述空间范围的上下边界值为一个长方体，所述空间范围为所述处理器计算域的rangebox；rangebox方案示意图如图5所示。

s107，根据所述处理器所有物质点邻域所覆盖的空间范围确定各个所述处理器相互重叠的范围；

s108，根据各个所述处理器相互重叠的范围确定重叠节点；

s109，假定k＝0，初始化tk时刻的物质点坐标；

s110，假定k＝0，采用局部最大熵插值函数计算tk时刻的物质点的局部形变数据；所述局部形变数据包括物质点形函数取值、物质点形函数导数取值、物质点体积、物质点密度、物质点变形梯度和物质点温度；

s111，假定k＝0，确定各个所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度；

s112，根据所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度确定全局质量矩阵、全局节点力和全局加速度；

本发明采用最优输运理论对力学平衡方程和能量平衡方程进行时间离散，该时间离散在模块s113到模块s114中实现，具体如下：

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

步骤2)：通过空间分散得到完全离散化的力学平衡方程和能量平衡方程，并采用分离算法进行求解；

步骤3)：假设温度场固定，显式求解力学平衡方程，得到材料的变形；

步骤4)：假设材料未有变形，采用牛顿法隐式求解能量平衡方程，得到材料的温度场。

采用了交错算法对强热流固耦合问题进行解耦计算，并采用了不同的算法求解动量守恒及能量守恒。其中，材料的动量守恒方程用显式算法求解，能量守恒方程采用隐式算法求解，保证了计算精度和计算效率。

s113，根据所述全局质量矩阵、全局节点力和全局加速度，采用最优输运理论对力学平衡方程进行时间离散，更新tk+1时刻的节点坐标。

s114，根据所述tk时刻的物质点的局部形变、所有节点在tk+1时刻的坐标，采用牛顿迭代法确定tk+1时刻所述节点的温度场；

s114具体包括：

令子迭代步n＝0并初始化节点温度、容差，计算等效热驱动力；

若所述节点温度大于所述容差，则确定温度增量与温度场，并继续下一步子迭代；

若所述节点温度不大于所述容差，则得到下一时间步的温度场。

本模块采用了分布式并行算法，提高了计算效率。本模块中分布式并行算法核心如下：

在计算材料的温度场时，采用newton-raphson迭代算法，首先根据节点总数目和处理器数目对全局矩阵进行划分，确定每个处理器所需处理的全局带状子矩阵范围；然后在各处理器上计算该处理器的局部矩阵，并将该局部矩阵根据处理器数目划分为局部带状子矩阵，并将相应的局部带状子矩阵传递给相应处理器；最后在各处理器上将接收的局部带状子矩阵组合为全局带状子矩阵，并输入分布式线性方程求解器进行求解，得到材料的温度场。

在s114之前还包括：

判断求解温度场的时间步长是否为动态响应时间步长的整数倍；

若所述求解温度场的时间步长为动态响应时间步长的整数倍，则进行所述根据所述tk时刻的物质点的局部形变、所有节点在tk+1时刻的坐标，采用牛顿迭代法确定tk+1时刻的所述物质点和所述节点的温度场的步骤；

若所述求解温度场的时间步长不为动态响应时间步长的整数倍，则到所述根据在tk+1时刻所有节点的坐标，利用局部最大熵插值函数，确定在tk+1时刻物质点的局部变形数据的步骤。

s115，根据tk+1时刻所述节点的坐标确定更新后该处理器的rangebox，并根据所述处理器的rangebox确定各个所述处理器相互重叠的范围；

s116，根据在tk+1时刻所有节点的坐标，更新物质点的坐标、体积和密度，并确定在tk+1时刻物质点的局部变形数据，所述局部形变数据包括物质点形函数取值、物质点形函数导数取值、物质点变形梯度和物质点温度；计算所得的局部最大熵插值形函数，在边界上满足克罗内克属性，在计算中可以直接添加位移边界条件。同时，局部最大熵插值函数值及其导数通过牛顿迭代法高效求解，求解效率与邻域中节点的个数无关，避免了求解一般无网格插值函数的高计算量，提高了计算效率。

s117，根据在tk+1时刻所述物质点的局部变形数据确定tk+1时刻的节点的局部形变数据；所述节点的形变数据包括局部节点力、局部节点动量矩阵以及局部节点质量矩阵；

s118，判断k+1是否等于所述计算步数；

s119，若k+1等于所述计算步数，则根据所述物质点和节点从t0时刻到tk+1时刻形变数据和坐标变化以及温度场变化确定强热流固耦合问题的解，完成动态响应分析；

s120，若k+1不等于所述计算步数，则令k＝k+1，根据该时刻所述物质点的坐标、体积以及密度确定更新后的物质点邻域所覆盖的空间范围，并返回s112进入下一次迭代。

作为一个具体的实施例，求解强热流固耦合问题的方法过程如下：

对于d维的连续介质问题域ω，将几何模型ω离散为一组物质点集{θp，k，p＝1，2，...，m；k＝0，1，...，n}和一组节点集{xa，k，a＝1，2，...，n；k＝0，1，...，n}。其中下标p和a代表物质点与节点，m和n为物质点与节点的个数，k代表时间步的索引号，第k和k+1个时间步分别用tk和tk+1进行表示。

第一步，定义总计算步数及处理器数量pⁱ，i＝1，...，q，在初始时刻k＝0，根据处理器数量q对该问题域中的物质点进行划分并将划分后的物质点集分别发送到各个处理器pⁱ中；

第二步，根据tk时刻各处理器pⁱ上各个物质点与节点的坐标以及物质点的邻域范围定义，获取各个物质点的邻域节点；

第三步，根据第二步处理器pⁱ上各物质点邻域nh(θp，k)确定该处理器上所有物质点邻域所覆盖的空间范围，并取该空间范围的上下边界值定义一个长方体作为该处理器计算域的rangebox；

第四步，根据第三步定义的各个处理器rangebox，确定各个rangebox之间相互重叠的部分定义为shadowbox，并将处于shadowbox的节点定义为shadownode，由此得到各处理器pⁱ上的shadownode，并将它们存储至共享数据二维交换表格ck，该表格的行号与列号分别代表需要进行数据交换的两个处理器的序号，表格中的数据为需要交换的节点数据；

第五步，在各个处理器pⁱ上初始化tk时刻物质点坐标θp，k、采用局部最大熵插值函数计算物质点形函数取值na，k(θp，k)、物质点形函数导数取值物质点体积vp，k、物质点密度ρp，k、物质点变形梯度

第六步，在各处理器pⁱ单独计算该处理器节点的局部质量对角矩阵局部节点力

和局部加速度并将该处理器上所有shadownode的相关数据发送至数据交换表ck中，并接收shadownode在其他各个处理器上的计算结果，接收到的节点信息与局部标号相同的节点计算结果进行求和，得到全局质量矩阵mⁱ、全局节点力fⁱ和全局加速度aⁱ，完成数据同步；

第七步，更新节点数据：利用上一步的结果进行数据更新，采用最优输运理论进行时间离散，对节点进行k→k+1时刻的瞬态分析，显式计算所有节点在k+1时刻的坐标其中xk+1＝{xa，k+1，a＝1，2，...，n}表示所有节点坐标的集合；

第八步，第六步到第七步已经完成对恒温假设下运动方程的求解得到材料的变形，接下来将转入能量方程的并行化求解从而得到材料的温度场。求解温度场时将在能量方程中引入第六步和第七步中材料变形所引起的能量耗散。本发明通过newton-raphson迭代法计算温度场。算法如下：

令子迭代步n＝0并初始化节点温度容差ε，计算等效热驱动力如果计算温度增量与温度场并使k＝k+1，继续下一步子迭代；如果则得到下一时间步的温度场

由于实际计算中求解动态响应的时间步长与求解温度场的时间步长不一致，我们假设求解温度场的时间步长为动态响应时间步长的n倍。在计算温度场之前，首先确定本时间步是否求解温度场，即k是否为n的整数倍。如果是，则进行第九步，激活温度场求解器。如果不是，则转到第十二步。进行物质点的更新。

第九步，本发明中的温度场算法需要针对全局的线性方程组进行求解。根据现有节点数量确定全局线性方程组的维度，并根据处理器数目q将全局矩阵划分为q个全局带状子矩阵，每个处理器负责组装其中之一；

第十步，各个处理器pⁱ根据分配得到的物质点子集，计算该处理器上的局部矩阵该局部矩阵仅仅包含该处理器上的节点相关信息。然后将该局部矩阵划分为q个局部带状子矩阵，并将相应的局部带状子矩阵分配给相应的处理器；

第十一步，各处理器pⁱ接收来自其他各个处理器的局部带状子矩阵，并组合成全局带状子矩阵，然后将该带状子矩阵输入分布式线性方程求解器进行求解。最终经过newton-raphson迭代法求得节点温度场ta，k+1；

第十二步，根据更新的节点坐标确定确定tk+1时刻各个处理器pⁱ的rangebox范围，并根据所述处理器的rangebox确定各个所述处理器相互重叠的范围，重新搜索shadownode并构建数据交换集ck+1；

第十三步，根据第七步得到的节点坐标xk+1，利用局部最大熵插值函数得到物质点在k+1时刻的局部变形数据，包括：物质点增量变形梯度

和物质点变形梯度

第十四步，根据第十三步得到的物质点在k+1时刻的物质点增量变形梯度和物质点变形梯度，更新k+1时刻的节点力fk+1，节点动量矩阵lk+1、节点质量矩阵mk+1；

第十五步，根据第十四步得到物质点从k→k+1时刻运动和变形数据，更新物质点的坐标θp，k+1＝xa，k+1na(θp，k)、体积vp，k+1＝det(fp，k→k+1)vp，k与密度其中mp为物质点质量；

第十六步，根据第十五步计算得到的物质点坐标θp，k+1、体积vp，k+1与密度ρp，k+1，通过邻域定义重新计算物质点邻域nh(θp，k+1)并采用局部最大熵插值函数重新计算邻域内节点的形函数值na(θp，k+1)及导数值

第十七步，判断时间步tk+1，若k+1＝n，代表已经计算至最后一步，此时退出计算，最终得到物质点和节点从t0→tn时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析；若k+1≠n，则转入第七步进行下一次迭代计算，直至时间步k+1＝n为止，最终得到物质点和节点从t0→tn时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析。

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

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

(2)初始化计算参数，在初始时刻k＝0时，利用第(1)步中离散得到的节点集和物质点集，初始化节点坐标xa，k、节点形函数na(θp，k)、节点形函数导数节点力fk、节点动量矩阵lk、节点质量矩阵mk、节点温度ta，k，初始化物质点坐标xp，k、物质点邻域nh(θp，k)、物质点体积vp，k、物质点密度ρp，k、物质点变形梯度fp，k；

(3)本发明提出了变分结构来描述材料的能量变化率，该公式可以表示为：

在公式(1)中，为材料变形，t为材料温度，z为与材料能量耗散有关的内部变量，k为材料动能，w为等效能量，ρ0为材料密度，b为体力，为外加牵引力，q为材料内部热源，t0为参考温度，为施加在材料上的热流。其中，等效能量可以表示为：

其中a为材料的helmholtz自由能，n为材料的熵，φ^*为材料粘性耗散，ψ^*为材料的塑性耗散，χ为材料的热传导耗散，为衡量热传导的参数。

对于k到k+1这一时间间隔，状态参数已知。本发明使用增量变分结构来描述这一时间间隔内材料变化的能量最小值。

则在k+1时刻的材料状态参数tk+1，zk+1可由上述变分结构的极值来获取。在这一变分结构中，增量等效能量可以写为

假设状态参数tk+1，zk+1相互独立，鉴于增量变分结构φn中并未包括相关参数zk+1的梯度项，因此在计算中该内部变量可以通过求解增量变分结构中相关项的极值来获得。因此增量等效能量又可以写为

其中

本发明采用最优输运理论对增量变分结构中的动能项进行时间离散。根据该理论，在时间间隔[tk，tk+1]内的增量动能最小值可以表示为该时间间隔内初始密度和最终密度的wasserstein距离：

对变分结构(1)中的其他项，采用后向欧拉近似方法进行时间离散，可以得到增量变分结构的时间离散形式：

对该时间离散形式的增量变分结构取极值，即可获得在时间间隔[tk，tk+1]内动量守恒和能量守恒的时间离散形式。因此，材料的动态响应和温度场可由该变分结构的极值条件获得，即：

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

(4)第(3)步已经获得了变分结构的时间离散形式，接下来要对该变分结构进行空间离散。本发明采用了标准ritz-galerkin方法进行空点离散，局部最大熵插值函数被用作构建位移场和温度场的插值函数。

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

同时对于需满足条件：

na(x)≥0，a∈[1，n](11)

∑na(x)＝1(12)

∑na(x)xa＝x(13)

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

其中，z(x，λ^*(x))为定义的分区函数，且

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

引入插值函数之后，物质点在tk+1时刻的位置、温度、变形梯度和温度梯度可以用以下公式表述：

其中为物质点θp，k的邻域内节点集合。在引入局部最大熵插值函数之后，对时间分散的增量变分结构取机制，可以得到完全离散化的力学平衡方程和能量平衡方程，如下所示：

其中，为节点xa在时刻tk+1从物质点所获取的质量，节点的加速度则采用中心差分获得：

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

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

其中nh(xa，k)为节点xa，k计算节点力的有效领域，θp，k为tk时刻物质点的坐标，为物质点tk时刻的总变形梯度，为形函数导数，wp，k为物质点体积。更新节点力后，根据公式(26)更新节点动量矩阵和公式(27)更新节点质量矩阵：

(7)上述步骤对力学平衡方程进行了显式求解，接下来要对能量平衡方程进行求解。在能量平衡方程中，将材料内部由于变形产生的热量耗散及热传导定义为等效内部热量，将材料内部热源产生的热量及材料外部传入内部的热量定义为等效外部热量，两者可以表述为

本发明通过newton-raphson迭代法计算温度场。算法如下：

令子迭代步n＝0并初始化节点温度容差ε，计算等效热驱动力如果再计算温度增量与温度场并使k＝k+1，继续下一步子迭代；如果则得到下一时间步的温度场

(8)根据第(5)步得到物质点从k→k+1时刻运动和变形数据，更新物质点的坐标θp，k+1＝xa，k+1na(θp，k)、体积vp，k+1＝det(fp，k→k+1)vp，k与密度

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

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

以上实例为phmm方法的具体实现过程，其分布式多进程并行化步骤如图3所示：

对于d维的连续介质问题域ω，将几何模型ω离散为一组物质点集{θp，k，p＝1，2，...，m；k＝0，1，...，n}和一组节点集{xa，k，a＝1，2，...，n；k＝0，1，...，n}。其中下标p和a代表物质点与节点，m和n为物质点与节点的个数，k代表时间步的索引号，第k和k+1个时间步分别用tk和tk+1进行表示，如图9所示。

(1)第一步，定义总计算步数n，定义分布式处理器数量pⁱ，i＝1，...，p，在初始时刻k＝0根据处理器数量p，对连续介质问题域中的物质点进行划分并将划分后的物质点集分别发送到各个处理器pⁱ中。本发明采用物质点和节点相结合的方式对所求问题域进行离散，通过上述方式获得了问题域离散化的节点集和物质点集，在本步骤中即是要对物质点集进行划分，划分成不同的子集，并将不同的子集分配到不同的处理器上进行分布式计算。本发明采用metis算法对物质点集划分，后续计算在物质点上进行，不涉及网格。

(2)第二步，根据tk时刻pⁱ处理器上各个物质点与节点的坐标以及物质点的邻域范围获取各个物质点的邻域节点，其中dp，k为物质点的动态邻域大小，物质点大小通过以下公式确定：

dp，k＝δx×h。

其中δx为人工系数，取值范围为2.0～5.0，h为物质点的大小。可以通过初试网格得到物质点的初始体积，随后网格被抛弃，物质点的体积通过在每个时间步求解质量守恒方程动态计算而得，由此计算出物质点的大小。

(3)第三步，根据第二步各处理器上物质点邻域nh(θp，k)确定该处理器pⁱ上所有物质点邻域所覆盖的空间范围，并以该空间范围的上下边界将其近似为一个长方体(rangebox)该处理器上的所有节点可以表示为具体为：对pⁱ处理器上的各物质点，计算和其中lk代表pⁱ处理器上物质点邻域的下边界，uk代表pⁱ处理器上物质点邻域的上边界，定义以lk和uk为对角顶点的长方形(二维情况)或者长方体(三维情况)为rangebox，在rangebox中包含了该处理器物质点邻域内的所有节点

(4)第四步，根据第三步确定的rangebox，将各个rangebox之间相互重叠的部分定义为shadowbox，并将处于shadowbox的节点定义为共享节点(shadownode)，得到pⁱ的shadownode，将他们存储至共享数据二维交换表格ck，该表格的行号与列号分别代表需要进行数据交换的两个处理器的序号，表格中的数据为需要交换的节点数据；

(5)第五步，初始化tk时刻物质点坐标θp，k、物质点形函数na，k(θp，k)、物质点形函数导数物质点体积vp，k、物质点密度ρp，k、物质点变形梯度fp，k；

(6)第六步，在各处理器pⁱ上计算局部质量矩阵局部节点力和局部加速度

(7)第七步，各个处理器pⁱ将所包含的共享节点的局部节点力局部质量矩阵局部加速度发送至数据交换表ck中，并接收共享节点在其他各个处理器计算结果，完成数据同步，将接收到的节点信息与局部标号相同的节点计算结果进行求和，得到全局质量矩阵全局节点力和全局加速度在tk→tk+1,pⁱ计算更新节点运动学信息之后，对rangebox中的节点进行判断，如果该节点属于shadownode则将数据信息发送到数据交换表ck,同时从ck接收来自其他处理器的计算结果，进行节点运动学数据的组装，ck也在每个时间步得到更新；

(8)第八步，更新节点数据，利用第七步的结果进行数据更新，得到tk+1时刻的节点坐标和节点速度

(9)第九步，第八步已经完成对恒温假设下运动方程的显式求解得到材料的变形，接下来将转入能量方程的并行化求解从而得到材料的温度场。求解温度场时将在能量方程中引入材料变形所引起的能量耗散。本发明通过newton-raphson迭代法计算温度场。该算法需要求解一个全局线性方程，首先根据全局线性方程的维度(所有节点的数目)及处理器数目确定每个处理器需要处理的矩阵行数。如图6所示，全局矩阵维度为6，处理器数目为2，因此每个处理器需要处理3行内容。

(10)第十步，首先令子迭代步n＝0并初始化节点温度容差ε，计算等效热驱动力如果各个处理器pⁱ根据第八步分配的得到的物质点子集计算该处理器上的局部矩阵将该局部矩阵划分为p个局部带状子矩阵并将相应的局部带状子矩阵分配给相应的处理器。如图6所示，在每个处理器上都将计算得出的局部矩阵划分为两个局部带状子矩阵，处理器1将该处理器上的局部子矩阵传递给处理器2，反之亦然。

(11)第十一步，各处理器pⁱ接收来自其他各个处理器的局部带状子矩阵并组合成全局带状子矩阵然后将该带状子矩阵输入分布式线性方程组求解器进行求解。最终经过迭代求得节点温度场ta，k+1。如图6所示，每个处理器将接收的全部局部子矩阵组合为全局子矩阵，然后输入线性方程求解器。通过newton-raphson迭代法得到温度场。

(12)第十二步，根据第八步更新的节点坐标更新各rangebox大小，搜索并更新rangebox内包含的节点集，

(13)第十三步，根据第十二步确定了tk+1时刻各处理器上的节点值，定义rangebox的边界坐标值lk+1和uk+1得到更新，因此rangebox在每个时间步的迭代计算中将会动态更新，使得共享节点集ck+1中的共享节点随着计算的进行将会被动态更新；

(14)第十四步，tk+1时刻，更新物质点的坐标θp，k+1＝xa，k+1na(θp，k)、体积vp，k+1＝det(fp，k→k+1)vp，k与密度其中mp为物质点质量，并重新计算物质点邻域nh(θp，k+1)，并更新邻域内节点的形函数值na(θp，k+1)及导数值

(15)第十五步，判断时间步tk+1，若k+1＝n，代表已经计算至最后一步，此时退出计算，最终得到物质点和节点从t0→tn时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析；若k+1≠n，则转入第七步进行下一次迭代计算，直至时间步k+1＝n为止，最终得到物质点和节点从t0→tn时间段内的物理信息和动力学数据，包括变形、应力、密度、位移、速度、加速度、温度，完成材料的动态响应分析。

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

(1)本发明采用物质点与节点相结合的方式对几何模型进行离散，计算过程不涉及网格，完全避免了有限元方法由于固定网格所带来的问题，同时由于插值在节点进行而积分在物质点进行，解决了绝大部分无网格方法中存在的拉应力不稳定性；

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

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

(4)本发明采用了交错算法对强热流固耦合问题进行解耦计算，并采用了不同的算法求解动量守恒及能量守恒。本发明采用显式算法求解动量守恒方程，采用隐式算法求解能量守恒方程，保证了计算精度和计算效率。

(5)本发明采用了分布式并行化算法，将问题域离散成若干物质点集，每个处理器负责一个物质点集的计算。在本发明中，动量守恒在各处理器上使用显式算法求解，并同步共享节点的相关信息。同时本发明采用隐式算法求解能量守恒，在各个处理器上求解局部矩阵并组合为全局矩阵，然后使用分布式线性方程求解器求解，从而获得材料的温度场。

图11为本发明所提供的一种求解强热流固耦合问题的并行无网格系统流程示意图，如图11所示，本发明所提供的一种求解强热流固耦合问题的并行无网格系统，包括：第一获取模块1101、物质点集和节点集确定模块1102、第二获取模块1103、划分模块1104、邻域节点确定模块1105、空间范围确定模块1106、相互重叠的范围确定模块1107、重叠节点确定模块1108、初始化模块1109、tk时刻的物质点的局部形变数据确定模块1110、节点局部信息确定模块1111、节点全局信息确定模块1112、tk+1时刻的节点坐标更新模块1113、温度场确定模块1114、tk+1时刻rangebox更新模块1115、tk+1时刻物质点的局部变形数据确定模块1116、tk+1时刻的节点的形变数据确定模块1117、第一判断模块1118、强热流固耦合问题求解完成模块1119和空间范围更新模块1120。

第一获取模块1101用于获取d维的连续介质问题域；

物质点集和节点集确定模块1102用于将所述d维的连续介质问题域进行离散，得到物质点集和节点集；所述物质点集包括多个物质点；所述节点集包括多个节点；所述物质点用于表示某个点设定范围内的物质的物理信息；所述物理信息包括质量以及密度；所述节点用于表示某个点的动力学的信息；所述动力学的信息包括速度以及加速度；

第二获取模块1103用于获取所述d维的连续介质问题域的计算步数以及处理器的个数；

划分模块1104用于根据所述处理器的个数对所述物质点集进行划分，并将划分后的物质点集分别发送到各处理器中；

邻域节点确定模块1105用于根据tk时刻每一所述处理器上各个物质点与所述节点集的坐标以及所述物质点的邻域范围确定每个所述物质点的邻域节点；

空间范围确定模块1106用于根据所述物质点的邻域节点确定所述处理器上所有物质点邻域所覆盖的空间范围；所述空间范围的上下边界值为一个长方体，所述空间范围为所述处理器计算域的rangebox；

相互重叠的范围确定模块1107用于根据所述处理器所有物质点邻域所覆盖的空间范围确定各个所述处理器相互重叠的范围；

重叠节点确定模块1108用于根据各个所述处理器相互重叠的范围确定重叠节点；

初始化模块1109用于初始化k＝0时tk时刻的物质点坐标；

tk时刻的物质点的局部形变数据确定模块1110用于采用局部最大熵插值函数计算k＝0时tk时刻的物质点的局部形变数据；所述局部形变数据包括物质点形函数取值、物质点形函数导数取值、物质点体积、物质点密度物质点变形梯度和物质点温度；

节点局部信息确定模块1111用于确定各个所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度；

节点全局信息确定模块1112用于根据所述处理器上所有节点以及所述重叠节点的局部质量对角矩阵、局部节点力和局部加速度确定全局质量矩阵、全局节点力和全局加速度；

tk+1时刻的节点坐标更新模块1113用于根据所述全局质量矩阵、全局节点力和全局加速度，采用最优输运理论进行时间离散，更新tk+1时刻的节点坐标；

温度场确定模块1114用于根据所述tk时刻的物质点的局部形变、所有节点在tk+1时刻的坐标，采用牛顿迭代法确定tk+1时刻的所述节点的温度场；

tk+1时刻rangebox更新模块1115用于根据tk+1时刻所述节点的坐标确定更新后该处理器的rangebox，并根据所述处理器的rangebox确定各个所述处理器相互重叠的范围；

tk+1时刻物质点的局部变形数据确定模块1116用于根据在tk+1时刻所有节点的坐标，更新tk+1时刻物质点的坐标、体积和密度，并确定在tk+1时刻物质点的局部变形数据，所述局部形变数据包括物质点形函数取值、物质点形函数导数取值、物质点变形梯度和物质点温度；

tk+1时刻的节点的形变数据确定模块1117用于根据在tk+1时刻所述物质点的局部变形数据确定tk+1时刻的节点的局部形变数据；所述节点的形变数据包括局部节点力、局部节点动量矩阵以及局部节点质量矩阵；

第一判断模块1118用于判断k+1是否等于所述计算步数；

强热流固耦合问题求解完成模块1119用于若k+1等于所述计算步数，则根据所述物质点和节点从t0时刻到tk+1时刻形变数据和坐标变化以及温度场变化确定强热流固耦合问题的解，完成动态响应分析；

空间范围更新模块1120用于若k+1不等于所述计算步数，则使得k＝k+1，并返回节点全局信息确定模块，进入下一个迭代。

所述物质点集和节点集确定模块1102具体包括：初始节点确定单元、初始物质点确定单元、节点确定单元和物质点确定单元。

初始节点确定单元用于根据所述d维的连续介质问题域，采用不同形状的单元进行网格划分，生成初始节点，得到离散的节点；不同形状的单元包括三角形单元和四面体单元；

初始物质点确定单元用于根据所述单元的几何中心确定所述物质点；

节点确定单元用于将所述动力学的信息存储在所述节点上；

物质点确定单元用于将所述物理信息存储在所述物质点上。

所述温度场确定模块1114具体包括：等效热驱动力确定单元、温度增量与温度场确定单元和下一时间步的温度场确定单元。

等效热驱动力确定单元用于令子迭代步n＝0并初始化节点温度、容差，计算等效热驱动力；

温度增量与温度场确定单元用于若所述节点温度大于所述容差，则确定温度增量与温度场，并继续下一步子迭代；

下一时间步的温度场确定单元用于若所述节点温度不大于所述容差，则得到下一时间步的温度场。

本发明所提供的一种求解强热流固耦合问题的并行无网格系统，还包括：第二判断模块、执行温度场确定模块和跳过温度场确定模块。

第二判断模块用于判断求解温度场的时间步长是否为动态响应时间步长的整数倍；

执行温度场确定模块用于若所述求解温度场的时间步长为动态响应时间步长的整数倍，则进行所述根据所述tk时刻的物质点的局部形变、所有节点在tk+1时刻的坐标，采用牛顿迭代法确定tk+1时刻的所述物质点和所述节点的温度场的步骤；

跳过温度场确定模块，用于若所述求解温度场的时间步长不为动态响应时间步长的整数倍，则到所述根据在tk+1时刻所有节点的坐标，利用局部最大熵插值函数，确定在tk+1时刻物质点的局部变形数据的步骤。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。