CN117951973B - 一种稳定的弹塑性实体光滑粒子动力学模拟方法 - Google Patents
一种稳定的弹塑性实体光滑粒子动力学模拟方法 Download PDFInfo
- Publication number
- CN117951973B CN117951973B CN202410348874.XA CN202410348874A CN117951973B CN 117951973 B CN117951973 B CN 117951973B CN 202410348874 A CN202410348874 A CN 202410348874A CN 117951973 B CN117951973 B CN 117951973B
- Authority
- CN
- China
- Prior art keywords
- sph
- stress
- particles
- particle
- representing
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 239000002245 particle Substances 0.000 title claims abstract description 169
- 238000000034 method Methods 0.000 title claims abstract description 72
- 238000004088 simulation Methods 0.000 title claims abstract description 28
- 238000004364 calculation method Methods 0.000 claims abstract description 36
- 239000011159 matrix material Substances 0.000 claims abstract description 13
- 238000012937 correction Methods 0.000 claims abstract description 10
- 239000000463 material Substances 0.000 claims description 39
- 230000008859 change Effects 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 7
- 230000001133 acceleration Effects 0.000 claims description 4
- 230000007704 transition Effects 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 101100379079 Emericella variicolor andA gene Proteins 0.000 claims description 2
- 239000007787 solid Substances 0.000 abstract description 8
- 238000004422 calculation algorithm Methods 0.000 abstract description 7
- 230000008878 coupling Effects 0.000 abstract description 5
- 238000010168 coupling process Methods 0.000 abstract description 5
- 238000005859 coupling reaction Methods 0.000 abstract description 5
- 229910052782 aluminium Inorganic materials 0.000 description 17
- XAGFODPZIPBFFR-UHFFFAOYSA-N aluminium Chemical compound [Al] XAGFODPZIPBFFR-UHFFFAOYSA-N 0.000 description 17
- 230000010355 oscillation Effects 0.000 description 7
- XEEYBQQBJWHFJM-UHFFFAOYSA-N Iron Chemical compound [Fe] XEEYBQQBJWHFJM-UHFFFAOYSA-N 0.000 description 6
- 229910052751 metal Inorganic materials 0.000 description 6
- 239000002184 metal Substances 0.000 description 6
- 229910000831 Steel Inorganic materials 0.000 description 3
- 238000006073 displacement reaction Methods 0.000 description 3
- 239000010959 steel Substances 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000006835 compression Effects 0.000 description 2
- 238000007906 compression Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000009863 impact test Methods 0.000 description 2
- 239000007769 metal material Substances 0.000 description 2
- 239000008188 pellet Substances 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 229910000838 Al alloy Inorganic materials 0.000 description 1
- 239000004411 aluminium Substances 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 230000005489 elastic deformation Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种稳定的弹塑性实体光滑粒子动力学模拟方法，属于耦合算法技术领域，具体步骤如下：将实体结构离散成设定数量的带有物理信息的SPH粒子；当空间位置处发生形变计算SPH粒子的核函数导数的修正矩阵，并对核函数导数进行一致性修正；计算空间位置的变形梯度和变形梯度率，并计算当前构型下的柯西应力；将柯西应力转化为第一Piola‑Kirchhoff应力并利用黎曼近似关系更新第一Piola‑Kirchhoff应力。采用上述一种稳定的弹塑性实体光滑粒子动力学模拟方法，将固体黎曼解与完全拉格朗日SPH结合，形成了一种新的自适应拉格朗日SPH格式，无需引入人工应力和人工粘性，消除了拉伸不稳定性，并成功解决了不同算例需要人工调节可调谐参数的问题。
Description
技术领域
本发明涉及耦合算法技术领域，尤其是涉及一种稳定的弹塑性实体光滑粒子动力学模拟方法。
背景技术
金属材料是工程和科学领域常用的结构材料，其在不同强度外载荷下表现出弹性变形、塑性流动甚至大变形断裂等特性。目前，对于结构受力变形行为的数值研究存在多种方法。其中，光滑粒子流体动力学方法（Smoothed Particle Hydrodynamics，简称SPH）是一种无网格方法。SPH方法在离散求解守恒方程时不再依赖于背景网格，因此能够避免传统网格方法在材料发生严重局部塑性变形时可能出现的数值不稳定性问题。因此，SPH方法非常适用于求解结构材料的弹塑性力学行为。
目前在固体动力学中，SPH方法有两种常用的格式。一种是更新拉格朗日SPH（ULSPH）格式，其中相邻粒子的识别和核函数的计算依赖于材料的当前构型。在每个时间步骤中，粒子都有可能进入或退出支持域。另一种是完全拉格朗日（TLSPH）格式，其相邻粒子的识别和核函数的计算依赖于材料的初始构型。ULSPH方法存在一个显著的问题，即拉伸不稳定性。目前最常用的解决方法是Monaghan和Gray等人开发的人工压力/应力方法。然而，这种方法需要频繁调整与人工应力相关的可调参数，给计算带来了不确定性。相比之下，ULSPH方法基于参考构型评估核函数，能够从本质上消除核函数和应力空间相互作用产生的奇异性，因此不会出现拉伸不稳定性问题。
尽管TLSPH在固体力学研究中取得了成功，但在处理激波问题时，引入额外的人工粘性力以确保稳定性仍然是一种常见的做法。在各种人工粘性方法中，挑战在于调整可调参数，存在显著的局限。与人工应力方法相比，黎曼（Riemann）算法可能是更有利的解决方案。Riemann求解器善于处理非线性冲击波动态界面问题，最早由Parshikov等人将Riemann算法引入到固体SPH，并证明了其在高速冲击和岩石力学等领域的应用可行性。然而，目前常用的耦合算法为Riemann-ULSPH，该格式无法克服拉伸不稳定问题。此外，在现有的Riemann-SPH中缺乏考虑金属弹塑性转变产生的分裂波对Riemann压力波系的影响。
发明内容
本发明的目的是提供一种稳定的弹塑性实体光滑粒子动力学模拟方法，以解决上述技术问题。
为实现上述目的，本发明提供了一种稳定的弹塑性实体光滑粒子动力学模拟方法，具体步骤如下：
步骤S1：将实体结构离散成设定数量的SPH粒子，且每个SPH粒子携带关于实体材料的物理信息；
步骤S2：当空间位置处发生形变后，计算SPH粒子/>的核函数导数的修正矩阵，并对核函数导数进行一致性修正；
步骤S3：计算空间位置的变形梯度和变形梯度率，根据变形梯度和变形梯度率计算得到当前构型下的柯西应力；
步骤S4：将柯西应力转化为第一Piola-Kirchhoff应力；
步骤S5：HLLC黎曼求解器计算SPH粒子的应力和速度的黎曼解，利用黎曼近似关系更新第一Piola-Kirchhoff应力。
优选的，在步骤S1中，物理信息包括质量、密度、内能、速度、体积以及应力。
优选的，在步骤S2中，
核函数导数的修正矩阵的计算公式如下：
;
其中，为修正矩形，/>为初始构型下SPH粒子i的修正矩阵，/>为SPH粒子的初始位置，/>，j表示SPH粒子i支持域内的第j个邻域SPH粒子，/>为SPH粒子i的初始位置，/>为SPH粒子j的初始位置，/>表示SPH粒子i与SPH粒子j的初始距离；/>表示核函数，核函数为Wendland函数，/>表示SPH粒子的体积，/>表示初始构型下的第j个邻域SPH粒子的体积，/>为梯度算子，/>为SPH粒子i的梯度算子，/>为向量外积符号。
优选的，在步骤S3中，
空间位置的变形梯度和变形梯度率的计算公式如下：
;
其中，表示变形梯度，/>表示SPH粒子i变形梯度，/>表示SPH粒子的空间位置，/>和/>分别为变形后的SPH粒子i和SPH粒子j的空间位置；
为变形梯度率，/>表示SPH粒子j的变形梯度率，/>表示SPH粒子的速度，。
优选的，在步骤S3中，计算得到当前构型下的柯西应力过程如下：
步骤S31：计算SPH粒子i的速度梯度；
速度梯度计算公式如下：
;
其中，为速度梯度，/>为SPH粒子i的速度梯度；
步骤S32：计算SPH粒子i的应变率张量和旋转率张量；
计算公式如下：
;
其中，T表示矩阵转置符，和/>分别表示应变率张量和旋转率张量，/>为SPH粒子i的应变率张量，/>为SPH粒子i的旋转率张量；
步骤S33：计算SPH粒子i的偏应力；
偏应力更新计算公式如下：
;
其中，G为材料剪切模量，I为单位矩阵，和/>分别为偏应力和偏应力变化率，/>表示SPH粒子i的偏应力，/>表示SPH粒子i的偏应力变化率，/>表示矩阵的迹；
步骤S34：计算SPH粒子i的压力；
SPH粒子i的压力计算公式如下：
;
其中，表示压力，/>表示SPH粒子i的压力，/>为材料初始密度，/>为材料密度，，/>，/>表示材料比内能，/>、/>以及/>均为材料常数；
步骤S35：计算SPH粒子i的柯西应力；
柯西应力的计算公式如下：
;
其中，为SPH粒子i的柯西应力，/>为SPH粒子的柯西应力。
优选的，在步骤S4中，柯西应力转化为第一Piola-Kirchhoff应力的公式如下：
;
其中，J表示变形梯度的行列式，为第一Piola-Kirchhoff应力。
优选的，在步骤S5中，计算SPH粒子i的加速度和内能变化率，具体公式如下：
;
其中，表示SPH粒子i的加速度，m表示SPH粒子的质量，/>为SPH粒子j的质量；为SPH粒子i通过黎曼解转化后的第一Piola-Kirchhoff应力，为SPH粒子j通过黎曼解转化后的第一Piola-Kirchhoff应力，/>和/>分别为SPH粒子i和SPH粒子j的初始密度；
表示内能变化率，/>为SPH粒子i的第一Piola-Kirchhoff应力，为SPH粒子对应的黎曼速度解。
优选的，和/>的计算公式如下：
;
;
;
其中，为SPH粒子j的柯西应力，/>为黎曼应力解；
和/>分别为通过黎曼应力解转化后的SPH粒子i和SPH粒子j的柯西应力；表示SPH粒子j的变形梯度，/>和/>分别表示SPH粒子i和SPH粒子j的变形梯度的行列式。
优选的，黎曼应力解的具体过程如下：
步骤S51：求解特征波速；
具体表达式如下：
;
;
其中，a表示热动力学声速，c为特征波速；
步骤S52：计算波后流场参数的预测值，计算公式如下：
;
;
;
其中，、/>以及/>分别为接触波波速、材料密度标量以及偏应力标量的预测值，/>；
和/>分别表示接触界面左右两侧粒子沿连线方向上的起止柯西应力标量，
和/>分别表示接触界面左右两侧粒子沿连线方向上的起止材料密度，/>表示接触界面左右两侧粒子沿连线方向材料密度标量；
和/>分别表示接触界面左右两侧粒子沿连线方向上的起止速度，/>为表示接触界面左右两侧粒子沿连线方向的速度标量；
为接触波波速，/>和/>分别表示接触界面左右两侧粒子沿连线方向上的起止接触波波速，/>表示接触界面左右两侧粒子沿连线方向的接触波波速标量；
为接触界面左右两侧粒子沿连线方向的偏应力标量；
步骤S53：判断界面两侧粒子是否存在弹塑性状态转变；
若满足，其中，/>为材料的屈服应力，说明波前材料处于弹性状态而波后材料状态进入了塑性，进一步确定进入塑性状态的波后流场参数：
;
;
;
;
;
其中，、/>、/>、/>以及/>分别为波后流场的偏应力标量、材料密度、压力、柯西应力标量以及速度标量的真实值；
系数，/>，/>，/>表示接触界面左右两侧粒子的比内能；
若不满足，说明界面两侧粒子没有弹塑性状态的转变，此时弹性波波后流场参数的真实值等于原始流场参数；
;
步骤S54：利用波后的流场参数的真实值重新计算接触波两侧的星形区域参量值；
;
;
;
其中，、/>以及/>分别为星形区域的接触波波速、材料密度标量以及偏应力标量，函数/>的表达式为：
;
其中，；
步骤S55：计算压力黎曼解；
计算公式如下：
;
;
其中，为SPH粒子对的压力黎曼解；
步骤S56：计算偏应力黎曼解；
计算公式如下：
;
其中，/>为SPH粒子对偏应力黎曼解；
步骤S57：计算总应力黎曼解；
计算公式如下：
；
其中，为粒子对中间应力黎曼解。
优选的，求解速度黎曼解公式如下：
;
;
;
其中，上标/>表示切线方向；/>表示速度黎曼解。
因此，本发明采用上述一种稳定的弹塑性实体光滑粒子动力学模拟方法，有益效果为：
通过将重新构建的弹塑性固体黎曼解与完全拉格朗日SPH（TLSPH）耦合，巧妙地结合了这两种方法的优势，创造出了一种有效而稳健的弹塑性Riemann-TLSPH耦合算法。该方法成功地克服了光滑粒子流体动力学（SPH）方法中固有的拉伸不稳定性和数值振荡问题。最重要的是，这种方法摆脱了对可调参数的人工应力方法的依赖，从而解决了其局限性。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法的完全拉格朗日SPH粒子接触示意图；
图2为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例1中的金属长杆自由端拉伸算例；
图3为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例1中的不同模拟方法金属长杆的应力场分布图；
图4为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例1中的不同模拟方法得到的金属长杆自由端速度剖面图；
图5为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例1中的不同模拟方法得到的金属长杆自由端位移剖面图；
图6为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例2中的铝梁动态冲击示意图；
图7为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例2中的冲击速度20m/s下铝梁不同时刻t的塑性应变云图；（a）t = 0.1 ms，（b）t = 0.5 ms，（c）t = 1.0ms，（d）t = 1.5 ms；
图8为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例2中的冲击速度20m/s下铝梁中段挠度时程曲线；
图9为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例2中的不同冲击速度下铝梁中段计算挠度与理论解对比图；
图10为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例3中的泰勒杆冲击示意图；
图11为本发明一种稳定的弹塑性实体光滑粒子动力学模拟方法实施例3中的Armco铁在t=60us不同视角下的变形模式；（a）2D视图，（b）3D视图。
附图标记
1、铝梁；11、固定端；2、钢制圆柱形弹丸。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
在本发明的描述中，需要说明的是，术语“上”、“下”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，或者是该发明产品使用时惯常摆放的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
在本发明的描述中，还需要说明的是，除非另有明确的规定和限定，术语“设置”、“安装”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
下面结合附图，对本发明的实施方式作详细说明。
实施例1
如图1所示，一种稳定的弹塑性实体光滑粒子动力学模拟方法的完全拉格朗日SPH粒子接触示意图，图中也示出了粒子对接触黎曼示意图。
本实施例以金属长杆自由端拉伸为算例，用于验证本发明的有效性。如图2所示，一个长度为L=10 m、横截面边长为H=0.2 m的铝合金杆。杆的左端固定，右端自由。杆的密度为2700kg/m³，弹性模量为72GPa，泊松比为0.3。状态方程采用Mie–Grüneisen方程，参数设定为C0 = 5328 m/s，= 2.0和S1＝1.338。初始条件包括杆的四分之一长度部分以初始速度v=5m/s向右移动，产生应力波。为了评估本实施例提出的方法的可靠性，采用了有限元方法（FEM）和引入人工粘性力的TLSPH方法作为对比，其中人工粘性力选用一阶vonNeumann和Richtmyer人工粘性，可调谐系数为/>。有限元分析使用非线性动力学软件LSDYNA进行，采用沙漏控制以减轻数值振荡，线性体积粘度系数和二次体积粘度系数设置为1.5。
图3展示了沿横截面边长粒子数为9，使用不同的计算方法，在t=0.06s时杆的应力分布。结果表明，有限元法计算得到的应力光滑性最差，应用人工粘性的TLSPH法次之，而本发明所提方法得到了高度平滑的应力分布。如图4-5所示，定量给出了铝杆自由端中心粒子的速度和位移时程曲线。有限元的速度剖面显示出显著的振荡，而施加人工粘性力的TLSPH方法同样存在明显的震荡。其结果与不同系数的选择有关。即使系数增加到1.5，仍然存在一些震荡。位移曲线显示了类似的结果。然而，与传统具有人工粘性力的TLSPH相比，本发明所提算法无需手动调整参数，并且实现了速度和位移剖面的完全无振荡。
实施例2
本实施例以铝梁1的动态冲击变形算例，验证本发明的有效性。如图6所示，铝梁1的几何尺寸为142.24毫米，厚度和宽度均为6.35毫米。铝梁1的两端为固定端11，一个长50毫米、直径14.74毫米的钢制圆柱形弹丸2以一定速度撞击铝梁1的中部。模拟中保持射弹与目标铝梁1的质量比和冲量比不变。将三维问题简化为二维平面问题，相关的计算材料参数见表2.1。
表2.1 铝梁受冲击变形算例计算参数
;
在图7中，展示了冲击速度为v=20 m/s时不同时刻铝梁1的塑性应变累积云图。从图7中可以观察到，局部塑性应变主要分布在两端固定位置和铝梁1的中段。整个冲击过程中，铝梁1未发生因拉伸不稳定性而导致的非物理断裂，验证了本发明所提方法能够克服拉伸不稳定性。图8展示了冲击速度为v=20 m/s时梁中心挠度的时程曲线。能够观察到梁的弹性恢复，可以预测在经过长时间的瞬态横向振荡后，铝梁最终达到永久塑性变形状态。为进一步定量验证，将计算的铝梁中心挠度与理论解析解进行了比较。在相同的工况设置下，对不同冲击速度进行了计算，其结果与理论解的比较见图9。模拟的挠度值与基于刚塑性分析的解析解吻合良好，验证了本发明所提方法的准确性。
实施例3
本实施例以泰勒杆冲击算例，验证验证本发明的有效性。
经典的泰勒杆冲击试验通常被用作数值方法的验证案例。如图10所示，长度为H、直径为d的钢制圆柱形弹丸2以一定速度撞击刚性壁，冲击压缩引发撞击体塑性变形。模拟了Johnson和Holmquist的泰勒杆冲击试验。圆柱形金属的长度和直径分别为25.46毫米和7.6毫米，材料为Armco铁，冲击速度为221 m/s。计算中采用Johnson–Cook本构模型的强度模型，具体材料参数见表3.1和表3.2。计算粒子间距为0.38 mm。同时为了模拟固体边界条件，在碰撞边界处应用镜像对称条件。
表3.1 Armco铁的材料常数
;/>
表3.2 Armco铁的Mie–Grüneisen状态方程参数
;
图11展示了60微秒时刻变形圆柱体的不同视图。模拟结果表明，Armco铁圆柱体的表面直径从上端到冲击端逐渐增加，主要变形集中在冲击端附近。这些变形特征与Johnson和Holmquist的实验观测结果相当一致。表3.3对模拟结果进行了定量比较，考虑到Armco铁圆柱体的情况，最终长度H1与初始长度H0的比值为0.737，最终直径d1与初始直径d0的比值为1.726。长度比和高度比的相对误差均在5%以内，与实验观测值吻合良好。因此，本发明提出的方法成功地再现了圆柱体的实验变形特性。
表3.3 不同方法得到的冲击体变形后的特征尺寸
;
综上所述，本发明采用了一种稳定的弹塑性实体光滑粒子动力学模拟方法，能够有效地模拟结构的动态拉伸和动态压缩变形。与传统方法相比，这种方法无需引入人工应力或人工粘性力，同时成功克服了传统SPH结构模拟中常见的拉伸不稳定性和数值不稳定性问题，确保了计算的准确性。这为金属材料动态变形模拟提供了一种可靠且高效的解决方案，为相关领域的研究和应用提供了有力支持。
最后应说明的是：以上实施例仅用以说明本发明的技术方案而非对其进行限制，尽管参照较佳实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对本发明的技术方案进行修改或者等同替换，而这些修改或者等同替换亦不能使修改后的技术方案脱离本发明技术方案的精神和范围。
Claims (5)
1.一种稳定的弹塑性实体光滑粒子动力学模拟方法，其特征在于，具体步骤如下：
步骤S1：将实体结构离散成设定数量的SPH粒子，且每个SPH粒子携带关于实体材料的物理信息；
步骤S2：当空间位置处发生形变后，计算SPH粒子/>的核函数导数的修正矩阵，并对核函数导数进行一致性修正；
在步骤S2中，
核函数导数的修正矩阵的计算公式如下：
；
其中，为修正矩形，/>为初始构型下SPH粒子i的修正矩阵，/>为SPH粒子的初始位置，/>，j表示SPH粒子i支持域内的第j个邻域SPH粒子，/>为SPH粒子i的初始位置，/>为SPH粒子j的初始位置，/>表示SPH粒子i与SPH粒子j的初始距离；/>表示核函数，核函数为Wendland函数，/>表示SPH粒子的体积，/>表示初始构型下的第j个邻域SPH粒子的体积，/>为梯度算子，/>为SPH粒子i的梯度算子，/>为向量外积符号；
步骤S3：计算空间位置的变形梯度和变形梯度率，根据变形梯度和变形梯度率计算得到当前构型下的柯西应力；
在步骤S3中，
空间位置的变形梯度和变形梯度率的计算公式如下：
；
；
其中，表示变形梯度，/>表示SPH粒子i变形梯度，/>表示SPH粒子的空间位置，/>和/>分别为变形后的SPH粒子i和SPH粒子j的空间位置；
为变形梯度率，/>表示SPH粒子j的变形梯度率，/>表示SPH粒子的速度；
在步骤S3中，计算得到当前构型下的柯西应力过程如下：
步骤S31：计算SPH粒子i的速度梯度；
速度梯度计算公式如下：
；
其中，为速度梯度，/>为SPH粒子i的速度梯度；
步骤S32：计算SPH粒子i的应变率张量和旋转率张量；
计算公式如下：
；
；
其中，T表示矩阵转置符，和/>分别表示应变率张量和旋转率张量，/>为SPH粒子i的应变率张量，/>为SPH粒子i的旋转率张量；
步骤S33：计算SPH粒子i的偏应力；
偏应力更新计算公式如下：
；
其中，G为材料剪切模量，I为单位矩阵，和/>分别为偏应力和偏应力变化率，/>表示SPH粒子i的偏应力，/>表示SPH粒子i的偏应力变化率，/>表示矩阵的迹；
步骤S34：计算SPH粒子i的压力；
SPH粒子i的压力计算公式如下：
；
其中，表示压力，/>表示SPH粒子i的压力，/>为材料初始密度，/>为材料密度，，/>，/>表示材料比内能，/>、 />以及/> 均为材料常数；
步骤S35：计算SPH粒子i的柯西应力；
柯西应力的计算公式如下：
；
其中，为SPH粒子i的柯西应力，/>为SPH粒子的柯西应力；
步骤S4：将柯西应力转化为第一Piola-Kirchhoff应力；
在步骤S4中，柯西应力转化为第一Piola-Kirchhoff应力的公式如下：
；
其中，J表示变形梯度的行列式，为第一Piola-Kirchhoff应力；
步骤S5：HLLC黎曼求解器计算SPH粒子的应力和速度的黎曼解，利用黎曼近似关系更新第一Piola-Kirchhoff应力；
在步骤S5中，计算SPH粒子i的加速度和内能变化率，具体公式如下：
；
；
其中，表示SPH粒子i的加速度，m表示SPH粒子的质量，/>为SPH粒子j的质量；为SPH粒子i通过黎曼解转化后的第一Piola-Kirchhoff应力，为SPH粒子j通过黎曼解转化后的第一Piola-Kirchhoff应力，/>和/>分别为SPH粒子i和SPH粒子j的初始密度；
表示内能变化率，/>为SPH粒子i的第一Piola-Kirchhoff应力，/>为SPH粒子对应的黎曼速度解。
2.根据权利要求1所述的一种稳定的弹塑性实体光滑粒子动力学模拟方法，其特征在于：在步骤S1中，物理信息包括质量、密度、内能、速度、体积以及应力。
3.根据权利要求2所述的一种稳定的弹塑性实体光滑粒子动力学模拟方法，其特征在于：
和/>的计算公式如下：
；
；
；
；
；
其中，为SPH粒子j的柯西应力，/>为黎曼应力解；
和/>分别为通过黎曼应力解转化后的SPH粒子i和SPH粒子j的柯西应力；/>表示SPH粒子j的变形梯度，/>和/>分别表示SPH粒子i和SPH粒子j的变形梯度的行列式。
4.根据权利要求3所述的一种稳定的弹塑性实体光滑粒子动力学模拟方法，其特征在于：黎曼应力解的具体过程如下：
步骤S51：求解特征波速；
具体表达式如下：
；
；
其中，a表示热动力学声速，c为特征波速；
步骤S52：计算波后流场参数的预测值，计算公式如下：
；
；
；
其中，、/>以及/>分别为接触波波速、材料密度标量以及偏应力标量的预测值，；
和/>分别表示接触界面左右两侧粒子沿连线方向上的起止柯西应力标量，
和/>分别表示接触界面左右两侧粒子沿连线方向上的起止材料密度，/>表示接触界面左右两侧粒子沿连线方向材料密度标量；
和/>分别表示接触界面左右两侧粒子沿连线方向上的起止速度，/>为表示接触界面左右两侧粒子沿连线方向的速度标量；
为接触波波速，/>和/>分别表示接触界面左右两侧粒子沿连线方向上的起止接触波波速，/>表示接触界面左右两侧粒子沿连线方向的接触波波速标量；
为接触界面左右两侧粒子沿连线方向的偏应力标量；
步骤S53：判断界面两侧粒子是否存在弹塑性状态转变；
若满足，其中，/>为材料的屈服应力，说明波前材料处于弹性状态而波后材料状态进入了塑性，进一步确定进入塑性状态的波后流场参数：
；
；
；
；
；
其中，、/>、/>、/>以及/>分别为波后流场的偏应力标量、材料密度、压力、柯西应力标量以及速度标量的真实值；
系数，/>，/>，/>表示接触界面左右两侧粒子的比内能；
若不满足，说明界面两侧粒子没有弹塑性状态的转变，此时弹性波波后流场参数的真实值等于原始流场参数；
；
步骤S54：利用波后的流场参数的真实值重新计算接触波两侧的星形区域参量值；
；
；
；
其中，、/>以及/>分别为星形区域的接触波波速、材料密度标量以及偏应力标量，函数/>的表达式为：
；
其中，
；
步骤S55：计算压力黎曼解；
计算公式如下：
；
；
其中，为SPH粒子对的压力黎曼解；
步骤S56：计算偏应力黎曼解；
计算公式如下：
；
其中，/>为SPH粒子对偏应力黎曼解；
步骤S57：计算总应力黎曼解；
计算公式如下：
；
其中，为粒子对中间应力黎曼解。
5.根据权利要求4所述的一种稳定的弹塑性实体光滑粒子动力学模拟方法，其特征在于，求解速度黎曼解公式如下：
；
；
；
；
其中，上标/>表示切线方向；/>表示速度黎曼解。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410348874.XA CN117951973B (zh) | 2024-03-26 | 2024-03-26 | 一种稳定的弹塑性实体光滑粒子动力学模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410348874.XA CN117951973B (zh) | 2024-03-26 | 2024-03-26 | 一种稳定的弹塑性实体光滑粒子动力学模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117951973A CN117951973A (zh) | 2024-04-30 |
| CN117951973B true CN117951973B (zh) | 2024-06-07 |
Family
ID=90802027
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410348874.XA Active CN117951973B (zh) | 2024-03-26 | 2024-03-26 | 一种稳定的弹塑性实体光滑粒子动力学模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117951973B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104200091A (zh) * | 2014-08-27 | 2014-12-10 | 北京航空航天大学 | 一种基于寿命模型分散性的柱塞泵加速因子区间计算方法 |
| CN111709197A (zh) * | 2020-06-17 | 2020-09-25 | 福州大学 | 基于黎曼不变量的sph入流边界处理方法 |
| CN112861445A (zh) * | 2020-12-29 | 2021-05-28 | 中国船舶重工集团公司第七研究院 | 一种无网格数值模型的模拟方法 |
| CN116992796A (zh) * | 2023-09-27 | 2023-11-03 | 中国科学技术大学 | 一种自适应低耗散的sph-hllc黎曼求解器耦合算法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2824598A1 (en) * | 2012-03-06 | 2015-01-14 | Fujitsu Limited | Simulation program, simulation method, and simulation device |
-
2024
- 2024-03-26 CN CN202410348874.XA patent/CN117951973B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104200091A (zh) * | 2014-08-27 | 2014-12-10 | 北京航空航天大学 | 一种基于寿命模型分散性的柱塞泵加速因子区间计算方法 |
| CN111709197A (zh) * | 2020-06-17 | 2020-09-25 | 福州大学 | 基于黎曼不变量的sph入流边界处理方法 |
| CN112861445A (zh) * | 2020-12-29 | 2021-05-28 | 中国船舶重工集团公司第七研究院 | 一种无网格数值模型的模拟方法 |
| CN116992796A (zh) * | 2023-09-27 | 2023-11-03 | 中国科学技术大学 | 一种自适应低耗散的sph-hllc黎曼求解器耦合算法 |
Non-Patent Citations (2)
| Title |
|---|
| Smoothed particle hydrodynamics simulation of converging Richtmyer-Meshkov instability;shenghong huang 等;《Physics of Fluids》;20200803;第32卷(第8期);第1-12页 * |
| 圆柱形汇聚激波诱导Richtmyer-Meshkov不稳定的SPH模拟;徐建于 等;《力学学报》;20190731;第51卷(第4期);第998-1011页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117951973A (zh) | 2024-04-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Geelen et al. | A phase-field formulation for dynamic cohesive fracture | |
| Borden et al. | A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects | |
| Johnson et al. | Another approach to a hybrid particle-finite element algorithm for high-velocity impact | |
| Ren et al. | Meshfree simulations of spall fracture | |
| Owen | ASPH modeling of material damage and failure | |
| Hollkamp | Modeling vibratory damage with reduced-order models and the generalized finite element method | |
| Dou et al. | A generalized plasticity model incorporating stress state, strain rate and temperature effects | |
| Zhou et al. | Coupled bending and torsional vibrations of non-uniform thin-walled beams by the transfer differential transform method and experiments | |
| Sobhani | Improvement of vibrational characteristics of multipurpose structures (plate and shells) used in aerospace components by deploying Graphene Oxide Powders (GOPs) in a matrix as a nano-reinforcement: A comprehensive study | |
| CN111488670B (zh) | 一种非线性的质点弹簧软组织形变仿真方法 | |
| Vasilev et al. | Experimental verification of a crystal plasticity-based simulation framework for predicting microstructure and geometric shape changes: Application to bending and Taylor impact testing of Zr | |
| CN117951973B (zh) | 一种稳定的弹塑性实体光滑粒子动力学模拟方法 | |
| Sadeghi et al. | Dynamic plastic behaviour of strain rate sensitive tubes under axial impact | |
| Li et al. | A 9-node co-rotational quadrilateral shell element | |
| Zhou et al. | A vector form conjugated-shear bond-based peridynamic model for crack initiation and propagation in linear elastic solids | |
| Patni et al. | Efficient modelling of beam-like structures with general non-prismatic, curved geometry | |
| Azzara et al. | Time response stress analysis of solid and reinforced thin-walled structures by component-wise models | |
| Hu et al. | Explicit phase-field material point method with the convected particle domain interpolation for impact/contact fracture in elastoplastic geomaterials | |
| Chen et al. | Accurate analytical approximation to post-buckling of column with Ramberg− Osgood constitutive law | |
| CN111639419B (zh) | 一种船舱内爆小变形塑性毁伤模式问题的判断方法 | |
| Yang et al. | Viscoelasticity dependence on hydrodynamic responses during water entry | |
| Tian et al. | Development of element model subroutines for implicit and explicit analysis considering large deformations | |
| Chen et al. | A coupled thermo-mechanical model for simulating the material failure evolution due to localized heating | |
| CN110334459A (zh) | 一种输电塔塔线体系快速精细化建模系统及方法 | |
| Chen et al. | Co-rotational formulations for geometrically nonlinear analysis of beam-columns including warping and Wagner effects |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |