TOTAL LAGRANGIAN MATERIAL POINT METHOD FOR THE DYNAMIC ANALYSIS OF NEARLY INCOMPRESSIBLE SOFT MATERIALS
-
摘要: 物质点法(MPM)在模拟非线性动力问题时具有很好的效果, 其已被广泛应用于许多大变形动力问题的分析中. 然而传统的MPM在模拟不可压或近似不可压材料的动力学行为时会产生体积自锁, 极大地影响模拟精度和收敛性. 本文针对近似不可压软材料的大变形动力学行为, 提出一种混合格式的显式完全拉格朗日物质点法(TLMPM). 首先基于近似不可压软材料的体积部分应变能密度, 引入关于静水压力的方程; 之后将该方程与动量方程基于显式物质点法框架进行离散, 并采用完全拉格朗日格式消除物质点跨网格产生的误差, 提升大变形问题的模拟精度; 对位移和压强场采用不同阶次的B样条插值函数并通过引入针对体积变形的重映射技术改进了算法, 提升算法的准确性. 此外, 算法通过实施一种交错求解格式在每个时间步对位移场和压强场依次进行求解. 最后, 给出几个典型数值算例来验证本文所提出的混合格式TLMPM的有效性和准确性, 计算结果表明该方法可以有效处理体积自锁, 准确地模拟近似不可压软材料的大变形动力学行为.
-
关键词:
- 完全拉格朗日物质点法 /
- 软材料 /
- 体积自锁 /
- 大变形 /
- 显式动力学
Abstract: The material point method (MPM) shows good performance in modeling nonlinear dynamic problems and has been widely used to simulate various types of large deformation dynamic problems. However, the classical MPM may suffer from the volumetric locking when modeling the dynamic responses of the incompressible or nearly incompressible materials, which reduces the computational accuracy and affects the convergence behavior greatly. In this work, a displacement-pressure mixed total Lagrangian material point method (TLMPM) with explicit time integration is proposed for the large deformation dynamic behavior of nearly incompressible soft materials. In this method, an equation about the hydrostatic pressure is introduced based on the volumetric part of the strain energy density of nearly incompressible soft materials. Then the introduced equation as well as the momentum equation is discretized within the framework of the explicit MPM and the total Lagrangian formulation is implemented to overcome the cell-crossing noise, which increases the computational accuracy for the problems involving large deformation. Furthermore, the B-spline interpolation functions with different orders are applied for displacement and pressure fields respectively and the mixed TLMPM is improved to increase the accuracy by introducing a remapping technique for the volumetric deformation. In addition, the staggered solving scheme is adopted and the displacement and the pressure are required to be solved sequentially in a single time step. Finally, several typical numerical examples are simulated by the mixed TLMPM and the convergence and accuracy are analyzed. The results demonstrate that the proposed mixed TLMPM is able to deal with the volumetric locking effectively and simulate the dynamic behavior of nearly incompressible soft materials involving large deformation accurately. -
引 言
软材料广泛存在于自然界和工程中, 如水凝胶、形状记忆聚合物等一些非线性弹性体, 以及人体或动物的组织和器官等一些生物材料[1]. 软材料往往具有高弹性、低模量、耐磨性、抗震性等优良特性, 并且在经历极大的变形时没有能量耗散[2], 因此在生活和工程中得到了广泛的应用. 对于这些软材料力学行为的分析主要包括理论[3]、实验[4]和数值模拟[5] 3种, 其中数值模拟相较理论方法适用范围更广, 并且不受实验环境和设备的制约[6]. 因此, 针对软材料力学行为分析的数值方法在近年来得到持续地发展和广泛地应用.
有限单元法FEM作为模拟固体力学行为最常用的一种数值方法, 被广泛地运用于模拟软材料的各种复杂变形[7-8]. 由于软材料往往具有不可压或近似不可压特性, 因此在处理该类材料时需要对算法进行特殊处理. 常见的处理方法包括位移-压强混合格式、B-bar和F-bar方法等[9]. 然而这些方法大多都用于静力问题的计算, 而对于不可压材料的显式动力学模拟, 目前相关的算法还比较少. 此外, 由于软材料在受外界刺激或载荷作用时可能会产生极端的大变形[10], 采用FEM进行模拟时可能会产生网格畸变, 导致模拟结果的误差较大或者难以收敛, 甚至可能模拟失败[11]. 因此, FEM在模拟大变形动力学问题时仍然具有一定的局限性.
为了更准确、有效地模拟大变形动力学问题, Sulsky等[12-13]于1994年将流体力学中的质点网格法加以改进, 提出了一种新的粒子类方法, 即物质点法(MPM). MPM将求解对象离散成物质点并设置独立于物质点的背景网格, 这种处理方式结合了拉格朗日和欧拉描述的优势, 使其在避免FEM中可能出现的网格畸变的同时也无需处理欧拉描述中的对流项. 因此近年来被广泛应用于高速冲击[14]、裂纹扩展[15]、边坡失稳[16]等诸多大变形强非线性问题的模拟中.
然而, 传统的MPM在模拟大变形问题时物质点会跨越背景网格, 其所采用的线性插值函数的梯度不连续性会产生较大的跨网格误差[17]. 为了提升MPM的精度, 众多学者对MPM进行了一系列改进, 提出了广义插值物质点法[18]、对流粒子域插值物质点法[19]和B样条物质点法[20]等方法. de Vaucorbeil等[21]基于B样条插值提出了一种完全拉格朗日物质点法TLMPM, 可以很好地消除这一误差并减小计算量. 该方法已经被成功应用于模拟大变形准静态[22]和动力冲击问题[23].
与FEM类似, MPM在计算不可压材料时同样会产生体积自锁, 也需要进行特殊处理. Kularathna等[24]与Zhang等[25]将Chorin提出的投影方法应用到MPM中, 用于求解不可压流体; Coombs等[26]将F-bar方法扩展到MPM中用于克服不可压材料中出现的体积自锁现象; Iaconeta等[27]提出了一种位移-压强混合格式的隐式MPM用于求解不可压材料的准静态大变形行为. 然而目前在求解不可压软材料的大变形动力学问题时, 显式MPM仍然缺乏有效的途径克服体积自锁.
本文在TLMPM理论框架下, 基于近似不可压软材料的体积应变能形式, 通过引入压强控制方程并进行离散, 提出一种位移-压强混合格式的显式TLMPM, 用于稳定求解近似不可压软材料的大变形动力学问题. 此外还给出了若干算例, 以验证所提出算法的准确性和有效性.
1. TLMPM基本理论
1.1 控制方程
对于大变形问题, 其初始构型上坐标为
$ \boldsymbol{X} $ 的点在变形后移动到位置$ \boldsymbol{x} $ . 可以据此计算变形梯度为$$ \boldsymbol{F} = \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{X}} = {\boldsymbol{1}} + {\nabla }_{0}\boldsymbol{u} $$ (1) 其中,
$ {\boldsymbol{1 }}$ 为2阶单位张量,$ \boldsymbol{u} = \boldsymbol{x}-\boldsymbol{X} $ 表示位移,$ \nabla $ 为梯度算子, 下标0表示相对于初始构型的量. 对于动力学问题, 其拉格朗日格式的动量方程为[28]$$ {\rho }_{0}\ddot{\boldsymbol{u}} = {\nabla }_{0}\cdot \boldsymbol{P} + {\boldsymbol{b}}_{0} $$ (2) 其中,
$ \rho $ 为密度,$ \ddot{\boldsymbol{u}} $ 表示$ \boldsymbol{u} $ 对时间的2阶导数,$ \boldsymbol{P} $ 为PK1应力,$ \boldsymbol{b} $ 为体力. 引入位移试函数$ {\boldsymbol{w}}^{u} $ 并应用高斯散度定理, 推导可得等效积分弱形式为$$ {\int }_{{\mathrm{\varOmega }}_{0}}{\rho }_{0}\ddot{\boldsymbol{u}}\cdot {\boldsymbol{w}}^{u}\mathrm{d}V = {\int }_{{\mathrm{\varGamma }}_{0}}{\stackrel-{\boldsymbol{t}}}_{0}\cdot {\boldsymbol{w}}^{u}\mathrm{d}S+ $$ $$\qquad{\int }_{{\mathrm{\varOmega }}_{0}}{\boldsymbol{b}}_{0}\cdot {\boldsymbol{w}}^{u}\mathrm{d}V-{\int }_{{\mathrm{\varOmega }}_{0}}\boldsymbol{P}:{\nabla }_{0}{\boldsymbol{w}}^{u}\mathrm{d}V $$ (3) 式中右端第1项为应力边界条件,
$ \stackrel{-}{\boldsymbol{t}} $ 表示面力,$ {\mathrm{\varOmega }}_{0} $ 和$ {\mathrm{\varGamma }}_{0} $ 分别表示积分域以及积分域的表面.1.2 离散及变量更新
与在计算过程中不断更新构型的传统MPM不同, TLMPM在求解力学问题时不参考当前构型, 而是将初始构型上的求解域离散成若干个物质点, 离散后物体的质量、体积等属性集中在物质点上. 同时在求解域范围内设置背景网格, 于是背景网格的等效质量可以由物质点得到, 即
$$ {m}_{I} = \sum _{p = 1}^{{n}_{p}}{m}_{p}{N}_{I}^{u}\left({\boldsymbol{X}}_{p}\right) = \sum _{p = 1}^{{n}_{p}}{\rho }_{0}^{p}{N}_{I}^{u}\left({\boldsymbol{X}}_{p}\right){V}_{0}^{p} $$ (4) 其中,
$ {n}_{p} $ 为物质点的总数, 上标和下标上的p均表示物质点p上的量, 下标$ I $ 表示背景网格节点上的量, 而$ {N}_{I}^{u}\left({\boldsymbol{X}}_{p}\right) $ 则表示初始构型上节点$ I $ 处的位移插值函数在物质点$ {\boldsymbol{X}}_{p} $ 处的取值. 而对于物质点上的试函数, 可以通过背景网格上的试函数进行插值映射获得, 即$$ {\boldsymbol{w}}_{p}^{u} = \sum _{I = 1}^{{n}_{I}}{N}_{I}^{u}\left({\boldsymbol{X}}_{p}\right){\boldsymbol{w}}_{I}^{u} $$ (5) $$ {\nabla }_{0}\cdot {\boldsymbol{w}}_{p}^{u} = \sum _{I = 1}^{{n}_{I}}{\boldsymbol{w}}_{I}^{u}\cdot {\nabla }_{0}{N}_{I}^{u}\left({\boldsymbol{X}}_{p}\right) $$ (6) 其中
$ {n}_{I} $ 为背景网格节点的总数. 类似地, 位移、速度和加速度也同样可以采用类似的方法进行映射.将插值映射格式带入弱形式方程(3), 经整理可得每个节点的离散控制方程为
$$ {m}_{I}{\ddot{\boldsymbol{u}}}_{I} = {\boldsymbol{F}}_{I}^{\mathrm{e}\mathrm{x}\mathrm{t}} + {\boldsymbol{F}}_{I}^{\mathrm{i}\mathrm{n}\mathrm{t}} $$ (7) 其中,
$ {\boldsymbol{F}}_{I}^{\mathrm{e}\mathrm{x}\mathrm{t}} $ 和$ {\boldsymbol{F}}_{I}^{\mathrm{i}\mathrm{n}\mathrm{t}} $ 分别为节点等效外力和内力, 它们可以通过下式计算$$ {\boldsymbol{F}}_{I}^{\mathrm{e}\mathrm{x}\mathrm{t}} = \sum _{p = 1}^{{n}_{p}}{N}_{I}^{u}\left({\boldsymbol{X}}_{p}\right){\boldsymbol{b}}_{0}^{p} + \sum _{p = 1}^{{n}_{p}}{N}_{I}^{u}\left({\boldsymbol{X}}_{p}\right){\stackrel-{\boldsymbol{t}}}_{0}^{p}{h}_{0}^{-1}{V}_{0}^{p} $$ (8) $$ {\boldsymbol{F}}_{I}^{\mathrm{i}\mathrm{n}\mathrm{t}} = -\sum _{p = 1}^{{n}_{p}}{\boldsymbol{P}}_{p}\cdot {\nabla }_{0}{N}_{I}^{u}\left({\boldsymbol{X}}_{p}\right){V}_{0}^{p} $$ (9) 其中
$ {h}_{0} $ 为边界层的厚度.在时间积分方面, 对于显式MPM, 可以采用向前差分和中心差分两种格式[21,29]. 为了简便计算, 考虑采用向前差分法进行积分, 便可得到时间、空间离散后的TLMPM离散求解格式, 即
$$ {m}_{I}\left({\dot{\boldsymbol{u}}}_{I}^{t + \Delta t}-{\dot{\boldsymbol{u}}}_{I}^{t}\right) = \Delta t\left({\boldsymbol{F}}_{I}^{\mathrm{e}\mathrm{x}\mathrm{t}} + {\boldsymbol{F}}_{I}^{\mathrm{i}\mathrm{n}\mathrm{t}}\right) $$ (10) 根据式(10)可以得到各个节点在
$ t + \Delta t $ 时刻的速度. 而对于物质点上的变量, 则需要根据节点上的变量进行逆映射, 其具体形式为$$ {\boldsymbol{u}}_{p}^{t + \Delta t} = {\boldsymbol{u}}_{p}^{t} + \Delta t\sum _{I = 1}^{{n}_{I}}{N}_{I}\left({\boldsymbol{X}}_{p}\right){\dot{\boldsymbol{u}}}_{I}^{t + \Delta t} $$ (11) $$ {\dot{\boldsymbol{u}}}_{p}^{t + \Delta t} = {\dot{\boldsymbol{u}}}_{p}^{t} + \Delta t\sum _{I = 1}^{{n}_{I}}{N}_{I}\left({\boldsymbol{X}}_{p}\right)\frac{{\boldsymbol{F}}_{I}^{\mathrm{e}\mathrm{x}\mathrm{t}} + {\boldsymbol{F}}_{I}^{\mathrm{i}\mathrm{n}\mathrm{t}}}{{M}_{I}} $$ (12) $$ {\boldsymbol{F}}_{p}^{t + \Delta t} = {\boldsymbol{1}} + \sum _{I = 1}^{{n}_{I}}{\nabla }_{0}{N}_{I}\left({\boldsymbol{X}}_{p}\right)\left({\boldsymbol{x}}_{I}^{t + \Delta t}-{\boldsymbol{X}}_{I}\right) $$ (13) 2. 位移-压强混合格式
2.1 压强控制方程及其离散
在对凝胶、生物材料等软材料进行数值求解时, 由于其往往具有不可压特性, 因此会产生体积自锁, 对计算结果的精度造成很大影响. 因此需要对当前的TLMPM进行改进, 以提升其精度.
在大变形框架下, 近似不可压软材料的应变能密度可以写成[30-31]
$$ W = {W}_{\mathrm{d}\mathrm{e}\mathrm{v}}\left(\boldsymbol{C}\right) + {W}_{\mathrm{v}\mathrm{o}\mathrm{l}}\left(J\right) $$ (14) 其中
$ \boldsymbol{C} = {\boldsymbol{F}}^{\mathrm{T}}\boldsymbol{F} $ ,$ {W}_{\mathrm{d}\mathrm{e}\mathrm{v}} $ 和$ {W}_{\mathrm{v}\mathrm{o}\mathrm{l}} $ 分别为应变能函数的偏量部分和体积部分. 这里可以假设软材料符合Neo-Hookean本构关系, 其偏量应变能密度$ {W}_{\mathrm{d}\mathrm{e}\mathrm{v}} $ 的形式为[9]$$ {W}_{\mathrm{d}\mathrm{e}\mathrm{v}}^{\mathrm{N}\mathrm{H}} = \frac{\mu }{2}\left({I}_{1}{I}_{3}^{-1/3}-3\right) $$ (15) 式中,
$ \mu $ 为剪切模量,$ {I}_{1} $ 和$I_3 $ 分别为为$ \boldsymbol{C} $ 的第一和第三不变量. 而对于近似不可压材料,$ {W}_{\mathrm{v}\mathrm{o}\mathrm{l}} $ 可以取$$ {W}_{\mathrm{v}\mathrm{o}\mathrm{l}} = \frac{K}{2}{\left(J-1\right)}^{2} $$ (16) 其中,
$ J = \mathrm{d}\mathrm{e}\mathrm{t}\left(\boldsymbol{F}\right) $ ,$ K $ 为体积模量.在采用TLMPM求解时为了避免体积自锁, 将应力分解成
$$ \boldsymbol{S} = {\boldsymbol{S}}_{\mathrm{d}\mathrm{e}\mathrm{v}} + pJ{\boldsymbol{C}}^{-1} = 2\frac{\partial {W}_{\mathrm{d}\mathrm{e}\mathrm{v}}}{\partial \boldsymbol{C}} + pJ{\boldsymbol{C}}^{-1} $$ (17) 其中
$ p $ 为静水压力. 它与体积应变能存在如下关系$$ \frac{\partial {W}_{\mathrm{v}\mathrm{o}\mathrm{l}}}{\partial J}-p = 0 $$ (18) 基于TLMPM框架将该方程离散求解, 便可对体应力进行修正以避免体积自锁, 从而得到改进后的混合TLMPM求解格式. 为了离散方程(18), 首先引入压强试函数
$ {w}^{q} $ 以推导出它的弱形式$$ {\int }_{{\mathrm{\varOmega }}_{0}}K\left(J-1\right){w}^{q}\mathrm{d}V = {\int }_{{\mathrm{\varOmega }}_{0}}p{w}^{q}\mathrm{d}V $$ (19) 与2.2节中位移试函数的离散类似, 物质点上的试函数
$ {w}_{p}^{q} $ 同样可以离散到背景网格上, 即$$ {w}_{p}^{q} = \sum _{I = 1}^{{n}_{I}}{N}_{I}^{q}\left({\boldsymbol{X}}_{p}\right){w}_{I}^{q} $$ (20) 式中,
$ {N}_{I}^{q}\left({\boldsymbol{X}}_{p}\right) $ 表示节点$ I $ 处的压强插值函数在物质点$ {\boldsymbol{X}}_{p} $ 处的取值. 将式(20)代入式(19), 得出离散格式的压强控制方程为$$ \left[\sum _{p = 1}^{{n}_{p}}{N}_{I}^{q}\left({\boldsymbol{X}}_{p}\right){V}_{p}^{0}\right]{p}_{I} = \sum _{p = 1}^{{n}_{p}}K\left({\stackrel-{J}}_{p}-1\right){N}_{I}^{q}\left({\boldsymbol{X}}_{p}\right){V}_{p}^{0} $$ (21) 式中
$ {\stackrel{-}{J}}_{p} $ 为等效体积比. 由于位移和压强插值需要采用的不同的插值函数, 因此$ {\stackrel{-}{J}}_{p} $ 不能直接使用物质点的变形梯度行列式, 而是需要将体积变形进行重映射. 首先需要将体积变形从物质点映射到节点上, 即$$ {\stackrel-{J}}_{I} = \frac{\displaystyle\sum _{p = 1}^{{n}_{p}}\mathrm{d}\mathrm{e}\mathrm{t}\left({\boldsymbol{F}}_{p}\right){N}_{I}^{q}\left({\boldsymbol{X}}_{p}\right){V}_{p}}{\displaystyle\sum _{p = 1}^{{n}_{p}}{N}_{I}^{q}\left({\boldsymbol{X}}_{p}\right){V}_{p}} $$ (22) 再将
$ {\stackrel{-}{J}}_{I} $ 从节点映射回物质点上, 得到$$ {\stackrel-{J}}_{p} = \sum _{I = 1}^{{n}_{I}}{N}_{I}^{q}\left({\boldsymbol{X}}_{p}\right){\stackrel-{J}}_{I} $$ (23) 于是, 物质点上的静水压力为
$$ {p}_{p} = \sum _{I = 1}^{{n}_{I}}{N}_{I}^{q}\left({\boldsymbol{X}}_{p}\right){p}_{I} $$ (24) 此外, 为了满足混合格式中的LBB条件[32]以保证算法的收敛性, 需要考虑采用不同阶次的B样条基函数分别对位移和压强进行插值. B样条函数可以较为简单地实现高阶插值, 并且易于改变插值的阶次, 用其作为插值函数可以提升MPM的精度, 实现混合格式的物质点法时无需增加新的节点, 程序实现较为简便[33]. 物质点法中具体的B样条插值实施细节可以参考文献[20]的相关工作. 在本文中, 如无其他说明, 位移和压强插值分别采用二次和线性B样条插值, 这种插值方法可以满足LBB条件, 保证算法的准确性和收敛性[34].
2.2 交错求解流程
为了同时考虑位移和压强场的求解, 混合格式的TLMPM采用交错求解的格式, 在一个时间步内依次进行位移和压强场的求解, 其流程如下.
(1) 位移场求解:
①将物质点的质量和体积映射到节点上;
②计算节点的内力、外力向量;
③求解位移场方程(7), 更新节点速度;
④更新物质点速度、位移、变形梯度;
(2) 压强场求解:
①将体积变形重映射, 计算物质点等效体积比;
②求解压强场方程(21), 更新节点压强;
③更新物质点的压强;
(3) 应用本构关系, 即式(17)更新应力.
3. 数值算例
为了验证本文提出的针对近似不可压软材料动力学的完全拉格朗日物质点算法, 本节给出了几个典型的数值算例. 在各个算例的模拟中, 为了消除时间步长对结果的影响, 时间步长统一设置为
$$ \Delta t = 0.1\mathrm{\Delta }{t}_{\mathrm{c}\mathrm{r}} $$ (25) 其中
$ \mathrm{\Delta }{t}_{\mathrm{c}\mathrm{r}} $ 为临界时间步长, 其确定方法可以参考文献[17].3.1 库克膜问题
首先考虑库克膜问题[34], 该问题通常被用来验证算法处理体积自锁的能力. 如图1所示, 模型左端固定, 右端受到大小为0.25 N/m的均匀剪力. 模型采用近似不可压材料, 其弹性模量和泊松比分别为1000 Pa和0.499, 密度为1 kg/m3. 为了验证算法的网格收敛性, 模拟时采用不同尺寸的背景网格并将模型离散成不同数量的物质点, 初始时每个背景网格内包含4个物质点, 采用标准格式和混合格式的TLMPM分别进行模拟, 计算获得的t = 3 s时模型右上角点的竖直位移如图2所示. 可以看出标准的TLMPM的模拟结果在物质点较密时逐渐趋近于混合格式的结果. 图3给出了两种算法在t = 3 s时的静水压力分布图, 从图中可以看出混合格式的 TLMPM相比标准格式能得到更光滑的结果. 这表明了本文提出的混合格式的TLMPM可以更有效地避免体积自锁, 得到更合理的应力结果.
![]()
图 2 库克膜: 采用不同数量物质点离散的t = 3 s时右上角竖向位移模拟结果Figure 2. Cook’s membrane: simulation results of the vertical displacement of the top right corner at t = 3 s using different numbers of material points![]()
图 3 库克膜: t = 3 s时的静水压力分布Figure 3. Cook’s membrane: hydrostatic pressure distribution at t = 3 s3.2 二维软梁的弯曲
考虑一个二维近似不可压软材料梁的大变形弯曲问题. 模型的尺寸、边界条件如图4所示, 梁的下端固定, 并整体给定初速度为10 m/s. 采用近似不可压Neo-Hookean材料, 对应的弹性模量和泊松比分别17 MPa和0.499, 密度为1100 kg/m3. 将模型采用不同密度的背景网格进行离散, 保持每个背景网格内包含4个物质点, 分别采用混合格式和标准格式的TLMPM进行模拟并考虑不同阶次的插值函数, 不同网格密度和插值函数的模拟结果分别如图5和图6所示. 图5(a)是不同网格尺寸的A点的位移时程曲线, 表明了混合格式的 TLMPM相比标准TLMPM收敛性更好. 对于网格尺寸为0.1 m的模型, 将插值函数的阶次降低, 分别得到线性位移插值的混合格式的 TLMPM和标准TLMPM模拟结果与之前的模拟进行对比, 其中线性位移插值的混合格式的 TLMPM为满足LBB条件压强采用0阶插值. 结果如图5(b)所示, 从图中可以看出混合格式的TLMPM即使降低了插值阶次仍然可以得到和二次位移插值较为相似的结果, 而线性插值的标准TLMPM结果与其他的相差较大. 表明了混合格式的 TLMPM在进行近似不可压软材料的大变形问题的模拟时比标准TLMPM具有更好的效果.
![]()
图 5 二维软梁的弯曲: A点x方向位移变化曲线Figure 5. Bending of 2D soft beam: history curves of x-displacement of point A![]()
图 6 二维软梁的弯曲: 不同时刻的静水压力分布. (a)~(b) 标准格式线性插值; (c)~(d) 标准格式二次插值; (e)~(f) 线性−常数混合格式; (g)~(h) 二次−线性混合Figure 6. Bending of 2D soft beam: hydrostatic distribution at different times. (a)~(b) Standard linear interpolation; (c)~(d) standard quadratic interpolation; (e)~(f) linear-constant mixed formulation and (g)~(h) quadratic-linear mixed formulation分别提取了4次模拟在t = 0.5, 1.0 s的静水压力分布, 如图6所示, 可以看出线性的标准TLMPM的构型与其他3种相差较大并且无法计算出合理的静水压力, 而混合格式 TLMPM在对近似不可压材料的体应力计算上可以得到理想的结果.
3.3 三维柱体扭转
最后, 考虑一个三维柱体的扭转问题. 一个方形截面柱底端受约束, 整个柱体受到与其坐标相关的初速度而发生扭转, 其构型如图7(a)所示, 柱体尺寸为1 m × 1 m × 6 m, 底面固定, 将坐标原点取在底面中心, 给柱体施加与位置坐标相关的初速度为:
${v}_{x} = -y\mathrm{sin}\left(\text{π} z/12\right)$ ,${v}_{y} = x\mathrm{sin}\left(\text{π} z/12\right)$ . 材料参数与3.2节的二维软梁相同. 将模型离散成20 × 20 × 120个物质点, 并设置网格尺寸为0.1 m, 即模型内部的每个背景网格内有8个物质点. 采用标准和混合格式的 TLMPM进行模拟, 得到0.1 s和0.3 s的静水压力分布如图7(b)~图7(e)所示. 从图中可以看出, 混合TLMPM在模拟三维大变形问题时, 静水压力结果在网格内和网格之间仍然只有很小的振荡, 而标准的TLMPM则受体积自锁影响, 无法得到合理的静水压力结果. 同时, 得到A点的z方向位移时程曲线如图8所示, 同时考虑了混合格式的传统MPM, 并将结果与改进单元技术的混合格式FEM结果[30]进行对比, 可以发现混合格式的TLMPM与FEM的结果相差很小, 而标准TLMPM和混合格式的MPM的结果与它们有一定差距.![]()
图 7 三维柱体扭转: 构型及静水压力分布. (a) 几何模型; (b) t = 0.1 s, 标准TLMPM; (c) t = 0.3 s, 标准TLMPM; (d) t = 0.1 s, 混合 TLMPM; (e) t = 0.3 s, 混合 TLMPMFigure 7. Twisting of a 3D column: configurations and distribution of hydrostatic pressure. (a) Geometry; (b) t = 0.1 s, standard TLMPM; (c) t = 0.3 s, standard TLMPM; (d) t = 0.1 s, mixed TLMPM and (e) t = 0.3 s, mixed TLMPM![]()
通过以上几个算例, 证明本文所提出的位移-压强混合格式的TLMPM算法在模拟近似不可压软材料的大变形动力学问题时可以有效防止体积自锁, 具有很好的收敛性和准确性.
4. 结论
本文基于MPM算法框架, 发展了一种位移-压强混合格式的完全拉格朗日物质点法(TLMPM)用于模拟近似不可压材料软材料的大变形动力学问题. 该方法通过体积应变能引入了关于静水压力的控制方程; 进一步基于MPM进行离散并采用完全拉格朗日格式以避免大变形问题中物质点跨网格所产生的误差; 对位移场和压强场进行交错求解, 并采用不同阶次的B样条基函数对位移和压强分别进行插值. 除此之外, 还针对体积变形引入重映射技术以保证静水压力求解的准确性. 通过几个二维和三维近似不可压问题的模拟, 验证了本文提出的混合格式的TLMPM可以有效避免体积自锁, 具有很好的收敛性并能够得到光滑的应力结果, 在模拟软材料大变形动力行为时效果相比标准TLMPM和混合格式的MPM都更为有效.
-
![]()
图 2 库克膜: 采用不同数量物质点离散的t = 3 s时右上角竖向位移模拟结果
Figure 2. Cook’s membrane: simulation results of the vertical displacement of the top right corner at t = 3 s using different numbers of material points
![]()
图 3 库克膜: t = 3 s时的静水压力分布
Figure 3. Cook’s membrane: hydrostatic pressure distribution at t = 3 s
![]()
图 5 二维软梁的弯曲: A点x方向位移变化曲线
Figure 5. Bending of 2D soft beam: history curves of x-displacement of point A
![]()
图 6 二维软梁的弯曲: 不同时刻的静水压力分布. (a)~(b) 标准格式线性插值; (c)~(d) 标准格式二次插值; (e)~(f) 线性−常数混合格式; (g)~(h) 二次−线性混合
Figure 6. Bending of 2D soft beam: hydrostatic distribution at different times. (a)~(b) Standard linear interpolation; (c)~(d) standard quadratic interpolation; (e)~(f) linear-constant mixed formulation and (g)~(h) quadratic-linear mixed formulation
![]()
图 7 三维柱体扭转: 构型及静水压力分布. (a) 几何模型; (b) t = 0.1 s, 标准TLMPM; (c) t = 0.3 s, 标准TLMPM; (d) t = 0.1 s, 混合 TLMPM; (e) t = 0.3 s, 混合 TLMPM
Figure 7. Twisting of a 3D column: configurations and distribution of hydrostatic pressure. (a) Geometry; (b) t = 0.1 s, standard TLMPM; (c) t = 0.3 s, standard TLMPM; (d) t = 0.1 s, mixed TLMPM and (e) t = 0.3 s, mixed TLMPM
-
[1] 彭向峰, 李录贤. 超弹性材料本构关系的最新研究进展. 力学学报, 2020, 52(5): 1221-1232 (Peng Xiangfeng, Li Luxian. State of the art of constitutive relations of hyperelastic materials. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(5): 1221-1232 (in Chinese) doi: 10.6052/0459-1879-20-189 [2] 郭辉, 胡文军, 陶俊林. 泡沫橡胶材料的超弹性本构模型. 计算力学学报, 2013, 30(4): 575-579 (Guo Hui, Hu Wen-jun, Tao Jun-lin. The superelasticty constitutive model for foam rubber materials. Chinese Journal of Computational Mechanics, 2013, 30(4): 575-579 (in Chinese) doi: 10.7511/jslx201304021 [3] 肖锐, 向玉海, 钟旦明等. 考虑缠结效应的超弹性本构模型. 力学学报, 2021, 53(4): 1028-1037 (Xiao Rui, Xiang Yuhai, Zhong Danming, et al. Hyperelastic model with entanglement effect. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(4): 1028-1037 (in Chinese) doi: 10.6052/0459-1879-21-008 [4] Liu R, Li Y, Liu Z. Experimental study of thermo-mechanical behavior of a thermosetting shape-memory polymer. Mechanics of Time-Dependent Materials, 2019, 23(3): 249-266
[5] 朱忠猛, 杨卓然, 蒋晗. 软材料黏接结构界面破坏研究综述. 力学学报, 2021, 53(7): 1807-1828 (Zhu Zhongmeng, Yang Zhuoran, Jiang Han. Review of interfacial debonding behavior of adhesive structures with soft materials. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(7): 1807-1828 (in Chinese) doi: 10.6052/0459-1879-21-131 [6] 冯西桥, 曹艳平, 李博. 软材料表面失稳力学. 北京: 科学出版社, 2017 Feng Xiqiao, Cao Yanping, Li Bo. Surface Instability Mechanics of Soft Materials. Beijing: Science Press, 2017 (in Chinese)
[7] 董金平, 张志强. 基于横观各向同性超弹性理论的短纤维增强橡胶本构模型的建立与应用. 计算力学学报, 2016, 33(2): 231-237 (Dong Jinping, Zhang Zhiqiang. Establishment and application of short fiber reinforced rubber constitutive model based on transversely isotropic hyperelastic theory. Chinese Journal of Computational Mechanics, 2016, 33(2): 231-237 (in Chinese) [8] 黄春阳, 唐山, 彭向和. 超弹性薄膜与可压缩基底双层结构表面失稳分析. 力学学报, 2017, 49(4): 758-762 (Huang Chunyang, Tang Shan, Peng Xianghe. Study of surface instability about hyperelastic films on auxetic substrates under compression. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(4): 758-762 (in Chinese) doi: 10.6052/0459-1879-17-161 [9] de Borst R, Crisfield MA, Remmers JC, et al. Non-Linear Finite Element Analysis of Solids and Structures. United Kingdom: John Wiley & Sons Ltd, 2012
[10] 郑晓静. 关于极端力学. 力学学报, 2019, 51(4): 1266-1272 (Zheng Xiaojing. Extreme mechanics. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(4): 1266-1272 (in Chinese) doi: 10.6052/0459-1879-19-189 [11] Hu ZQ, Liu Y, Zhang HW, et al. Implicit material point method with convected particle domain interpolation for consolidation and dynamic analysis of saturated porous media with massive deformation. International Journal of Applied Mechanics, 2021, 13(2): 2150023 doi: 10.1142/S175882512150023X
[12] Sulsky D, Chen Z, Schreyer HL. A particle method for history-dependent materials. Computer Methods in Applied Mechanics and Engineering, 1994, 118(1-2): 179-196 doi: 10.1016/0045-7825(94)90112-0
[13] Sulsky D, Zhou SJ, Schreyer HL. Application of a particle-in-cell method to solid mechanics. Computer Physics Communications, 1995, 87(1-2): 236-252 doi: 10.1016/0010-4655(94)00170-7
[14] Pan XF, Xu AG, Zhang GC, et al. Generalized interpolation material point approach to high melting explosive with cavities under shock. Journal of Physics D-Applied Physics, 2008, 41(1): 015401 doi: 10.1088/0022-3727/41/1/015401
[15] Zhang ZJ, Qiu YS, Hu ZQ, et al. Explicit phase-field total Lagrangian material point method for the dynamic fracture of hyperelastic materials. Computer Methods in Applied Mechanics and Engineering, 2022, 398: 115234 doi: 10.1016/j.cma.2022.115234
[16] 孙玉进. 岩土大变形问题的物质点法研究. 北京: 清华大学, 2017 (Sun Yujin. Research on Geotechnical Problems Involving Extremely Large Deformation Using The Material Point Method. Beijing: Tsinghua University, 2017 (in Chinese))
[17] 张雄, 廉艳平, 刘岩等. 物质点法. 北京: 清华大学出版社, 2013 Zhang Xiong, Lian Yanping, Liu Yan, et al. Material Point Method. Beijing: Tsinghua University Press, 2013 (in Chinese)
[18] Bardenhagen SG, Kober EM. The generalized interpolation material point method. Computer Modeling in Engineering & Sciences, 2004, 5(6): 477-495
[19] Sadeghirad A, Brannon RM, Burghardt J. A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations. International Journal for Numerical Methods in Engineering, 2011, 86(12): 1435-1456 doi: 10.1002/nme.3110
[20] Gan Y, Sun Z, Chen Z, et al. Enhancement of the material point method using B-spline basis functions. International Journal for Numerical Methods in Engineering, 2018, 113(3): 411-431
[21] de Vaucorbeil A, Nguyen VP, Hutchinson CR. A total-Lagrangian material point method for solid mechanics problems involving large deformations. Computer Methods in Applied Mechanics and Engineering, 2020, 360: 112783 doi: 10.1016/j.cma.2019.112783
[22] Zhang ZJ, Pan YX, Wang JH, et al. A total-Lagrangian material point method for coupled growth and massive deformation of incompressible soft materials. International Journal for Numerical Methods in Engineering, 2021, 122(21): 6180-6202 doi: 10.1002/nme.6787
[23] de Vaucorbeil A, Nguyen VP. Modelling contacts with a total Lagrangian material point method. Computer Methods in Applied Mechanics and Engineering, 2021, 373: 113503
[24] Kularathna S, Soga K. Implicit formulation of material point method for analysis of incompressible materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 673-686 doi: 10.1016/j.cma.2016.10.013
[25] Zhang F, Zhang X, Sze KY, et al. Incompressible material point method for free surface flow. Journal of Computational Physics, 2017, 330: 92-110 doi: 10.1016/j.jcp.2016.10.064
[26] Coombs WM, Charlton TJ, Cortis M, et al. Overcoming volumetric locking in material point methods. Computer Methods in Applied Mechanics and Engineering, 2018, 333: 1-21
[27] Iaconeta I, Larese A, Rossi R, et al. A stabilized mixed implicit Material Point Method for non-linear incompressible solid mechanics. Computational Mechanics, 2019, 63(6): 1243-1260 doi: 10.1007/s00466-018-1647-9
[28] 李锡夔, 郭旭, 段庆林. 连续介质力学引论. 北京: 科学出版社, 2015 Li Xikui, Guo Xu, Duan Qinglin. Introduction to Continuum Mechanics. Beijing: Science Press, 2015 (in Chinese)
[29] Kakouris EG, Triantafyllou SP. Phase-field material point method for dynamic brittle fracture with isotropic and anisotropic surface energy. Computer Methods in Applied Mechanics and Engineering, 2019, 357: 112503 doi: 10.1016/j.cma.2019.06.014
[30] Kadapa C. Novel quadratic Bézier triangular and tetrahedral elements using existing mesh generators: Extension to nearly incompressible implicit and explicit elastodynamics in finite strains. International Journal for Numerical Methods in Engineering, 2019, 119(2): 75-104 doi: 10.1002/nme.6042
[31] Tian FC, Zeng J, Zhang MN, et al. Mixed displacement–pressure-phase field framework for finite strain fracture of nearly incompressible hyperelastic materials. Computer Methods in Applied Mechanics and Engineering, 2022, 394: 114933 doi: 10.1016/j.cma.2022.114933
[32] Hughes TJR. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Englewood Cliffs, NJ: Prentice-Hall, 1987
[33] 王建华. 软物质溶胀行为分析的实体-壳元有限元法和实体元等几何法研究. [博士论文]. 大连: 大连理工大学, 2021 Wang Jianhua. Study on solid-shell finite element method and solid element isogeometric analysis method for the swelling deformation of soft materials. [PhD Thesis]. Dalian: Dalian University of Technology, 2021 (in Chinese)
[34] Telikicherla RM, Moutsanidis G. Treatment of near-incompressibility and volumetric locking in higher order material point methods. Computer Methods in Applied Mechanics and Engineering, 2022, 395: 114985 doi: 10.1016/j.cma.2022.114985
下载:
计量
- 文章访问数:
- HTML全文浏览量:
- PDF下载量:
