CN110457785B - 一种用于结构大变形响应的物质点法的物质信息映射方法 - Google Patents

Abstract

本发明提供一种用于结构大变形响应的物质点法的物质信息映射方法，该方法通过追踪结构大变形过程中物质点体积域的变化，建立物质点体积域与背景网格之间的拓扑关系；由背景网格节点将物质点的体积域划分为若干体积子域，在物质点各体积子域中心分别布置虚拟映射物质点，其中虚拟映射物质点上的物质信息由该物质点的物质点上的物质信息直接确定；进而基于虚拟映射物质点实现背景网格节点的质量、动量、体力、面力和应力等物质信息的映射。本发明能够准确给出物质点法可有效消除物质点法在求解结构大变形问题时的网格穿越误差，解决了物质点法在求解结构大变形动力学响应时求解精度不足的问题。

Description

一种用于结构大变形响应的物质点法的物质信息映射方法

技术领域

本发明涉及结构动力学技术领域，特别涉及一种用于结构大变形响应的物质点法的物质信息映射方法。

背景技术

物质点法是一种粒子型无网格算法，由一系列的拉格朗日型物质点离散求解域、采用欧拉背景网格覆盖整个求解域；其中，物质点携带材料的所有物质信息，包括位置、速度、动量、应力应变等历史变量，并跟随物体运动而运动；采用有限元法中的线性插值形函数实现各物质点和背景网格节点的信息交换，进而在背景网格上进行空间导数和控制方程的求解，实现相邻物质点之间的相互作用与联系。物质点法结合了拉格朗日法和欧拉法各自的优点，同时避免了各自的不足，具有显著的算法优势，已在各类结构大变形问题中得到了广泛应用。

然而，物质点法在求解背景网格节点内力时采用物质点积分，即物质点既作为材料离散点跟随结构运动而运动，又作为积分点以用于更新背景网格节点的内力信息。同时，物质点法采用线性插值形函数来实现物质点和背景网格节点之间的信息交换；通过线性插值形函数的一阶导数实现背景网格节点内力的求解。因此，当结构有较大变形时，物质点将不可避免地穿越背景网格边界，由于在背景网格边界处插值形函数的一阶导数具有突变，因而将引起背景网格节点内力的突变，而造成网格穿越误差，降低了物质点法的求解精度。

现有技术：或者通过提高物质点法的插值形函数阶次以减小物质点法的网格穿越误差。然而，随着插值形函数阶次的提高，不仅求解算法复杂，而且扩大了插值形函数的影响范围，导致在边界处难以满足插值特性即不满足Kronekerδ函数性质，带来了处理本质边界条件的困难。或者通过在背景网格内布置积分点以减小物质点法的网格穿越误差。但是，由于布置在背景网格内积分点的应力信息无法直接确定，因此需要额外的插值算法进行求解，增加了求解算法的复杂性。或者通过在物质点体积域内布置高斯积分点以用于求解背景网格节点的内力信息，但由于未考虑高斯积分点质量和动量等物质信息对背景网格节点的影响，造成背景网格节点质量和动量的不守恒，影响了求解精度。

发明内容

针对现有技术的不足，本发明提供一种用于结构大变形响应的物质点法的物质信息映射方法，该方法实现算法简单，不仅能够准确给出背景网格节点的内力信息，同时可满足背景网格节点动量和质量守恒，可有效消除物质点法在求解结构大变形问题时的网格穿越误差，提高了物质点法分析结构大变形动力学响应的求解精度。

为实现上述目的，本发明采取的技术方案如下：

一种用于结构大变形响应的物质点法的物质信息映射方法，包括如下步骤：

步骤1、基于拉格朗日物质点离散问题域，并初始化各物质点的质量、动量、密度、应力应变、变形梯度等物质信息；同时，布置一套规则的欧拉型背景网格以实现控制方程的求解和各物质点之间的相互作用和联系；

步骤2、通过物质点的变形梯度，确定各物质点的当前体积域；

步骤3、确定各物质点体积域与背景网格之间的拓扑关系，通过背景网格节点将物质点体积域划分为k个体积子域，在各个体积子域中心布置虚拟映射物质点，虚拟映射物质点上的物质信息由物质点上的物质信息求解得到；

步骤4、基于虚拟映射物质点，采用对应的线性插值形函数或其梯度，将虚拟映射物质点上的质量、动量、面力、体力和应力等物质信息映射到相应的背景网格节点上；

步骤5、在背景网格节点上施加本质边界条件，对控制方程进行显式求解；

步骤6、更新虚拟映射物质点的加速度、速度、应变增量、旋度增量和变形梯度等物质信息；

步骤7、基于虚拟映射物质点上加速度、速度、应变增量、旋度增量和变形梯度等物质信息，求解各物质点上相应的物质信息；

步骤8、更新物质点的位移、密度、体积和应力应变等物质信息；

步骤9、输出结果信息，进入下一计算流程或结束计算。

本发明解决了传统物质点法在求解结构大变形响应时的网格穿越误差，避免了由于插值形函数阶次的提高而引起的求解算法复杂和本质边界条件难以施加的问题，克服了由于在背景网格内布置积分点而带来的复杂插值算法的不足，同时避免了仅通过高斯积分点对背景网格节点内力进行积分而造成的背景网格节点质量和动量不守恒的问题，能够给出结构在动力学作用下更为精确的动态响应，为研究结构大变形动力学响应提供了一种高精度的分析方法。

附图说明

图1为本发明步骤3物质点体积域与背景网格之间的拓扑关系以及虚拟映射物质点的布置示意图；

图2为本发明步骤4线性插值形函数及其导数示意图；

图3为本发明实施例2的计算模型；

图4为本发明实施例2应力及位移沿轴线的分布及求解误差；

图5为本发明实施例3的计算模型；

图6为本发明实施例3的应力分布云图。

具体实施方式

实施例1

一种用于结构大变形响应的物质点法的物质信息映射方法，包括如下步骤：

步骤1、基于拉格朗日物质点离散问题域，并初始化各物质点的质量、动量、密度、应力应变、变形梯度等物质信息；同时，布置一套规则的欧拉型背景网格以实现控制方程的求解和各物质点之间的相互作用和联系；

步骤2、通过物质点的变形梯度，确定各物质点的当前体积域；

在t时刻，物质点p的体积域，可表示为：

其中：

为t时刻物质点p的变形梯度；

和

分别表示物质点p在初始时刻和t时刻的体积域。

步骤3、确定各物质点体积域与背景网格之间的拓扑关系，通过背景网格节点将物质点体积域划分为若干体积子域，在各个体积子域中心布置虚拟映射物质点，虚拟映射物质点上的物质信息由物质点上的物质信息求解得到；

首先，根据t时刻下物质点p的体积域

确定物质点p的体积域在x_i方向(i＝1,2,3)的分布范围

和

进而，通过各背景网格节点的坐标，确定物质点p的体积域与背景网格相交而形成的k个体积子域；最后，在物质点p的各个体积子域中心布置单个虚拟映射物质点p′，其位置信息记为x_p′，如图2所示。对N_p个物质点共布置得到N_p′个虚拟映射物质点，其中虚拟映射物质点p′所在体积子域的体积记为Ω_p′，其对应物质点p的体积记为V_p。则虚拟映射物质点p′的质量m_p′、动量

变形梯度F_p′，面力

体力b_p′和应力σ_p′等物质信息可由相应物质点p上的对应物质信息求得，分别为：

b_p′＝b_p；σ_p′＝σ_p

其中：m_p为物质点p的质量，在计算过程中保持不变即满足质量守恒；V_p为物质点p的体积；

b_p和σ_p分别为物质点p的面力、体力和应力信息。

步骤4、基于虚拟映射物质点，采用对应的线性插值形函数或其梯度，将虚拟映射物质点上的质量、动量、面力、体力和应力等物质信息映射到相应的背景网格节点上；

式中：

为t时刻第I个背景网格节点的质量；

为t时刻虚拟映射物质点p′对应在背景网格节点I上的插值形函数；

为虚拟物质点p′对应在背景网格节点I上的插值形函数梯度；

和

分别表示t时刻背景网格节点I和虚拟映射物质点p′的动量；

和

分别为t时刻背景网格节点I上的外力和内力；

为t时刻，作用在物体上的表面力，h为面力作用面的厚度；

为t时刻虚拟物质点p′的体积；

为t时刻虚拟物质点p′所受的体力。

步骤5、在背景网格节点上施加本质边界条件，对控制方程进行显式求解；

对于结构动力学问题，位移边界条件可以统一写为速度和加速度的形式，因此t时刻的本质边界条件为：

式中：

和

为t时刻本质边界

处背景网格节点I上的加速度；

和

为t时刻本质边界

处的背景网格节点I上的速度。

t时刻，背景网格节点I加速度的求解公式为：

t+Δt时刻，背景网格节点I上速度的显式积分格式为：

式中：Δt为显式求解时间步长。

步骤6、更新虚拟映射物质点的加速度

速度

应变增量Δε_p′、旋度增量ΔΩ_p′和变形梯度F_p′等物质信息，分别为：

步骤7、基于虚拟映射物质点上加速度、速度、应变增量、旋度增量和变形梯度等物质信息，求解各物质点上相应的物质信息；

其中，在步骤3中将物质点p的体积域划分为了k个体积子域，即物质点p对应与k个虚拟映射物质点p′，因此物质点p的物质信息可由k个虚拟映射物质点上的相应的物质信息计算得到：

式中：

和

分别为t时刻虚拟物质点p′和物质点p的体积；

Δε_p、ΔΩ_p和F_p分别为物质点p的速度、应变增量、旋度增量和变形梯度。

步骤8、更新物质点的位移、密度、体积和应力应变等物质信息；

t+Δt时刻物质点的位移u_p、密度ρ_p、体积V_p和应力σ_p、应变ε_p分别为：

式中：

表示材料的本构模型。

步骤9、输出结果信息，进入下一计算流程或结束计算。

根据需要输出各物质点的位置、动量、密度和应力应变等物质信息，进入下一计算流程，即重复步骤2-步骤9；或结束计算。

实施例2

本实施例以一维杆为例，具体说明本发明方法的应用方式；

如图3所示，具有初速度的、两端固定的一维杆，由于初速度的作用杆将产生沿轴向的动态响应。分别利用物质点法和本发明对此一维杆振动响应进行分析，计算参数选取统一为：弹性模量E＝50Pa，密度ρ＝25kg/m³，模型长度L＝1m，v₀＝0.1m/s，网格尺寸为6.25mm，每个背景网格内布置4个物质点，时间步长为0.5μs。图4分别给出了t＝0.5s时刻，物质点法和本发明所给出的位移和应力沿杆轴线的分布情况以及相应的L₂误差。

实施例3

本实施例以二维自由梁为例，具体说明本发明方法的应用方式；

如图5所示，长度L＝0.06m、高度H＝0.01m的等截面直梁，其中a点为梁中心点，梁的弹性模量为E＝318GPa、密度ρ＝1845kg/m3、泊松比v＝0.054。初始时刻t＝0，梁内各点分布有一阶自由梁弯曲振动模态速度；在t>0时，梁在自由状态下发生振动，分别利用物质点法和本发明对此二维自由梁的振动响应进行模拟，其中背景网格尺寸为dx＝dy＝6.25×10^-4m，时间步长dt＝0.001μs，每个网格内布置4×4个物质点。图6给出了t＝8μs时刻，物质点法和本发明计算所得梁在x方向的正应力云图。

Claims (5)

1.一种用于结构大变形响应的物质点法的物质信息映射方法，其特征在于，包括如下步骤：

步骤1、基于拉格朗日物质点离散问题域，初始化各物质点的物质信息，包括质量、动量、密度、应力应变、变形梯度；布置规则的欧拉型背景网格；

步骤2、通过物质点的变形梯度，确定各物质点的当前体积域；

步骤3、确定各物质点体积域与背景网格之间的拓扑关系，通过背景网格节点将物质点体积域划分为k个体积子域，在各个体积子域中心布置虚拟映射物质点，虚拟映射物质点上的物质信息由物质点上的物质信息求解得到，包括：

根据t时刻下物质点p的体积域

确定物质点p的体积域在x_i方向(i＝1,2,3)的分布范围

和

通过各背景网格节点的坐标，确定物质点p的体积域与背景网格相交而形成的k个体积子域；在物质点p的各个体积子域中心布置单个虚拟映射物质点p′，其位置信息记为x_p′；对N_p个物质点共布置得到N_p′个虚拟映射物质点，其中虚拟映射物质点p′所在体积子域的体积记为Ω_p′，其对应物质点p的体积记为V_p；则虚拟映射物质点p′的质量m_p′、动量

变形梯度F_p′，面力

体力b_p′和应力σ_p′由相应物质点p上的对应物质信息求得，分别为：

b_p′＝b_p；σ_p′＝σ_p

其中，m_p在计算过程中保持不变；

步骤4、基于虚拟映射物质点，采用对应的线性插值形函数或其梯度，将虚拟映射物质点上的物质信息，包括质量、动量、面力、体力和应力映射到相应的背景网格节点上；

映射到背景网格节点上的质量、动量、面力、体力和应力物质信息表示如下：

式中：

为t时刻第I个背景网格节点的质量；

为t时刻虚拟映射物质点p′对应在背景网格节点I上的插值形函数；

为虚拟物质点p′对应在背景网格节点I上的插值形函数梯度；

和

分别表示t时刻背景网格节点I和虚拟映射物质点p′的动量；

和

分别为t时刻背景网格节点I上的外力和内力；

为t时刻，作用在物体上的表面力，h为面力作用面的厚度；

为t时刻虚拟物质点p′的体积；

为t时刻虚拟物质点p′所受的体力；

步骤5、在背景网格节点上施加本质边界条件，对控制方程进行显式求解；

步骤6、更新虚拟映射物质点的加速度、速度、应变增量、旋度增量和变形梯度物质信息；

步骤7、基于虚拟映射物质点上加速度、速度、应变增量、旋度增量和变形梯度物质信息，求解各物质点上相应的物质信息；

步骤8、更新物质点的位移、密度、体积和应力应变物质信息；

步骤9：输出结果信息，重复步骤2-9或结束计算。

2.根据权利要求1所述的方法，其特征在于，所述步骤3中，体积子域划分方式为：根据t时刻下物质点的体积域，确定物质点的体积域在x_i方向(i＝1,2,3)的分布范围；进而通过各背景网格节点的坐标，确定物质点的体积域与背景网格相交而形成的体积子域。

3.根据权利要求2所述的方法，其特征在于，所述步骤3中，虚拟映射物质点上的物质信息由物质点上相应的物质信息、物质点的体积以及虚拟映射物质点所在体积子域的体积直接确定；背景网格节点的质量、动量和内力物质信息由虚拟映射物质点上的质量、动量、面力、体力和应力物质信息映射得到。

4.根据权利要求1所述的方法，其特征在于，所述步骤5中，位移边界条件写为速度和加速度的形式，则t时刻的本质边界条件为：

式中：

和

为t时刻本质边界

处背景网格节点I上的加速度；

和

为t时刻本质边界

处的背景网格节点I上的速度。

5.根据权利要求1所述的方法，其特征在于，所述步骤7中，物质点对应于k个虚拟映射物质点，物质点的物质信息由k个虚拟映射物质点上的相应的物质信息计算得到。

Images