一维理想弹塑性体的SPH-HLLC耦合算法

王展铭; 陈龙奎; 黄生洪

doi:10.11883/bzycj-2024-0004

一维理想弹塑性体的SPH-HLLC耦合算法

doi: 10.11883/bzycj-2024-0004

中国科学技术大学工程科学学院材料力学行为与设计重点实验室，安徽合肥 230026

详细信息

作者简介:
王展铭（1997－　），男，博士研究生，wangzmrr@mail.ustc.edu.cn

通讯作者:
黄生洪（1974－　），男，博士，研究员，博士生导师，hshnpu@ustc.edu.cn

中图分类号: O389
计量
- 文章访问数: 302
- HTML全文浏览量: 93
- PDF下载量: 90
- 被引次数: 0
出版历程
- 收稿日期: 2024-01-02
- 修回日期: 2024-05-14
- 网络出版日期: 2024-05-15
- 刊出日期: 2024-08-05

SPH-HLLC coupled method for one-dimentional elastic-perfectly plastic model

Key Laboratory of Mechanical Behavior and Design of Materials, School of Engineering Science, University of Science and Technology of China, Hefei 230026, Anhui, China

摘要

摘要: 通过弹塑性波分析求得HLLC（Harten-Lax-van Leer-contact）近似黎曼解，提出了SPH（smoothed particle hydrodynamics）与一维理想弹塑性体模型下近似的HLLC黎曼求解器耦合的一种构造简单的算法。在SPH计算中，支持域内每个粒子对都存在一个黎曼间断问题，它的黎曼解被代入控制方程中计算。其中一维理想弹塑性体的HLLC近似黎曼解的思想是：先假设整体处于弹性状态计算黎曼解，然后对计算结果进行塑性条件修正，最后用修正后的物理变量计算HLLC近似黎曼解。将提出的SPH-HLLC耦合算法与传统SPH算法在一维算例下的计算结果进行对比，结果表明，该算法能有效模拟一维理想弹塑性体材料的碰撞，并能有效抑制在不同材料之间的压强和偏应力震荡，这是传统SPH方法很难做到的。
- 弹塑性数值模拟 /
- 固体碰撞 /
- SPH黎曼耦合算法
Abstract: A 1D SPH (smoothed particle hydrodynamics) and approximate HLLC (Harten-Lax-van Leer-contact) Riemann solver coupled method for elastic-perfectly plastic model is proposed through elastic and plastic wave analysis. In SPH simulations, each particle pair in the supporting domain generates a Riemann problem, whose solutions are substituted into governing equations. The philosophy of HLLC approximate Riemann solver is to divide the procedure into three steps: assume the whole state in elastic deformation and compute Riemann problem, and then reconstruct flux under von Mises yielding conditions and compute the final HLLC Riemann solution with reconstructed fluxes. We compare the new SPH-HLLC method with the traditional SPH method in several numerical tests, which show that this method can effectively simulate collision and reflected rarefaction waves between the materials, and it can profoundly suppress oscillations of pressure and deviatoric stress at contact interface between different materials, which the traditional SPH method finds difficult to realize. Moreover, the new SPH-HLLC scheme shows better energy performance than the traditional SPH method in 2D test case where initial kinetic energy is successfully transformed into internal energy with new SPH-HLLC scheme while total energy significantly decreases with time using the traditional SPH method.
- elastic-plastic numerical simulation /
- solid collision /
- SPH-Riemann scheme

HTML全文

光滑粒子流体动力学（smoothed particle hydrodynamics，SPH）方法是一种构造简单的具有很好鲁棒性的无网格粒子方法，它将整个计算系统离散化为一系列带有自身物理属性的粒子，粒子的位移、压强、速度、密度等属性随时间推进，被广泛应用于各种工况场景计算之中，并取得了良好的计算结果，比如自由表面流动、多相界面流动、爆炸冲击现象、激波模拟、磁流体力学等。

为了维持计算稳定性，早期SPH在计算中需要添加人工黏性项^[1]，然而，这种方法需要根据不同工况手动调节人工黏性项的参数，为了维持稳定，有些还需要额外增加一些稳定项，具有一定的局限性。Vila^[2]、Parshikov等^[3]将黎曼求解器引入SPH框架中，较好地克服了上述人工黏性方法的缺陷，从此，SPH黎曼求解器耦合算法就被广泛应用于各种流体相关算例中，在多相流体运动中展现出超越传统SPH方法对于界面模拟的优势。然而，虽然后续人们对SPH黎曼求解器耦合算法进行了很多的研究，但是对于固体相互作用的SPH黎曼求解器耦合算法的研究，大多还是将黎曼求解器的液体解形式代入计算^[4-5]。这样的方法可以起到一定的替代人工黏性项的作用，但是对于波系的捕捉不太准确，因此，对于物质界面的处理也会出现一定的非物理震荡。

在网格算法中，已经存在一些包含弹塑性波分析的完全黎曼解。Lin等^[6]对理想弹塑性固体提出了一种基于迭代方法的近似黎曼算子。姚成宝等^[7]针对具有高度非线性特征通用形式的Mie-Grüneisen状态方程和流体弹塑性本构方程的多种介质间的相互作用问题，提出了一种通用、健壮的多介质黎曼问题求解方法。然而，这类迭代求得完全黎曼解的方式在SPH计算过程中比较繁琐耗时，我们希望就像在SPH模拟液体流动过程一样，在SPH框架下建立一个构造简单、计算方便的近似黎曼解，然而，SPH与包含波系分析的完全黎曼解或者近似黎曼解耦合的研究还没有出现。事实上，已经有很多学者对理想弹塑性固体的HLLC(Harten-Lax-van Leer-contact) 近似黎曼解进行了一定的研究^[8-11]，这为其代入SPH粒子对求解提供了必要的参考。在这些研究中，近似求解HLLC近似黎曼解被分为三步，首先，假设界面两侧左右状态尚未达到弹性极限，以左右皆为弹性波的三波系假设计算临时HLLC近似黎曼解；然后，根据已有的临时近似黎曼解进行Mises应力检测，如果至少有一侧产生塑性变形，则根据密度积分和Rankine-Hugoniot条件重构出塑性状态物理参数；最后，再对重构后的左右状态计算出最后的HLLC近似黎曼解。我们的SPH内理想弹塑性固体粒子对的HLLC近似黎曼解也遵循这一套逻辑求解。类似的波系分析逻辑可以扩展到二维三维流固耦合算例中形成统一黎曼解，而在传统算法中，这样的算例需要额外的镜像粒子防止流体粒子侵入固体界面。

本文中将提出一种理想弹塑性固体相互作用的一维SPH-HLLC耦合算法，在追求构造形式简单、计算方便的同时可以较为准确地捕捉弹塑性波，通过对比固体强冲击强激波工况下一维算例的理论解以及二维圆柱内碰撞算例，验证算法正确性。

1. 计算方法

1.1 SPH控制方程及其离散化

选用三次b样条核函数^[12]：

$W\left(q,h\right)={\alpha }_{v}\left\{ \begin{array}{ll}1-\dfrac{3}{2}{q}^{2}+\dfrac{3}{4}{q}^{3}& \qquad 0{\text{≤}} q{\text{≤}} 1\\ \dfrac{1}{4}{\left(2-q\right)}^{3}& \qquad 1{\text{≤}} q{\text{≤}} 2\\ 0& \qquad 其他\end{array}\right.$

(1)

在一维问题中，系数 ${\alpha _v} = \dfrac{1}{{6h}}$ 。 $q = {{| {{x_i} - {x_j}} |}/h}$ ，表示位于x_i和x_j的两粒子之间的归一化距离，其中h为光滑长度。

对于一个忽略热传导过程的系统来说，SPH离散后的控制方程如下：

$\left\{ \begin{gathered} \frac{{{\text{d}}{\rho _i}}}{{{\text{d}}t}} = \sum\limits_{j = 1}^n {{m_j}\left( {{\boldsymbol{u}}_i^\beta - {\boldsymbol{u}}_j^\beta } \right){\boldsymbol{L}}_i^{\beta \gamma }\frac{{\partial {W_{ij}}}}{{\partial {\boldsymbol{x}}_i^\gamma }}} \\[-1pt] \frac{{{\text{d}}{\boldsymbol{u}}_i^\alpha }}{{{\text{d}}t}} = \frac{2}{{{\rho _i}}}\sum\limits_{j = 1}^n {{\boldsymbol{\sigma }}_*^{\alpha \beta }{\boldsymbol{L}}_i^{\beta \gamma }\frac{{\partial {W_{ij}}}}{{\partial {\boldsymbol{x}}_i^\gamma }}} {V_j} \\[-1pt] \frac{{{\text{d}}{e_i}}}{{{\text{d}}t}} = \frac{2}{{{\rho _i}}}\sum\limits_{j = 1}^n {\left( {{\boldsymbol{u}}_*^\alpha - {\boldsymbol{u}}_i^\alpha } \right){\boldsymbol{\sigma }}_*^{\alpha \beta }{\boldsymbol{L}}_i^{\beta \gamma }\frac{{\partial {W_{ij}}}}{{\partial {\boldsymbol{x}}_i^\gamma }}{V_j}} \\[-1pt] {\boldsymbol{L}}_i^{\beta \gamma } = - {\left( {\sum\limits_{j = 1}^n {{\boldsymbol{x}}_{ij}^\beta \frac{{\partial {W_{ij}}}}{{\partial {\boldsymbol{x}}_i^\gamma }}} {V_j}} \right)^{ - 1}} \\ \end{gathered} \right.$

(2)

式中： $\rho$ 为SPH粒子的密度， $m$ 为质量， ${\boldsymbol{u}}$ 为速度矢量， ${\boldsymbol{L}}$ 为梯度修正项，W为SPH核函数， ${\boldsymbol{x}}$ 为粒子的笛卡尔坐标， ${\boldsymbol{\sigma }}$ 为应力张量， $V$ 为粒子体积， $e$ 为粒子的内能，下标i表示粒子编号，下标j表示粒子i支持域内的第j个粒子， ${{\boldsymbol{\sigma }}_*}$ 为应力的HLLC黎曼解，上标 $\alpha 、\beta$ 表示矢量或者张量的分量。而应力由球对称压力 $p$ 和偏应力张量 ${\boldsymbol{S}}$ 两部分组成， ${\boldsymbol{I}}$ 为单位矩阵，表示如下：

${\boldsymbol{\sigma }} = - p{\boldsymbol{I}} + {\boldsymbol{S}}$

(3)

一维条件下偏应力 $S$ 由如下方式更新：

$\left\{ \begin{gathered} \frac{{{\text{d}}{S_i}}}{{{\text{d}}t}} = \frac{{4\mu }}{3}\frac{{\partial {u_i}}}{{\partial {x_i}}} \\[-1pt] \left| S \right| {\text{≤}} \frac{2}{3}{Y_0} \\ \end{gathered} \right.$

(4)

式中： $\mu$ 为物质的剪切模量， ${Y_0}$ 为屈服极限。

1.2 一维SPH-HLLC黎曼解

为更简洁地研究黎曼求解器内的波系分化，进行如下的坐标转化^[13]，在拉格朗日坐标系下研究黎曼解问题：

${\text{d}}t = {\text{d}}\bar t,\qquad\quad {\text{d}}x = u{\text{d}}\bar t + \frac{1}{\rho }{\text{d}}\bar \xi$

(5)

式中： $\bar t$ 和 $\bar \xi$ 分别为转化后拉格朗日坐标系下的时间和物质坐标。在拉格朗日坐标系下，弹性波（上标为E）和塑性波（上标为P）的波速^[11]分别为：

${c^{\text{E}}} = \sqrt {\frac{{\partial p}}{{\partial \rho }} + \frac{p}{{{\rho ^2}}}\frac{{\partial p}}{{\partial e}} + \frac{{4\mu }}{\rho }} , \qquad\quad{c^{\text{P}}} = \sqrt {\frac{{\partial p}}{{\partial \rho }} + \frac{p}{{{\rho ^2}}}\frac{{\partial p}}{{\partial e}} - \frac{{3\mu }}{{\rho Y_0^2}}{S^2}}$

(6)

简化HLLC黎曼解分三步求解。首先，假设黎曼求解器左右状态只产生三个波：左行弹性波、接触波和右行弹性波。通过假设接触波两侧应力相等 ${\sigma _{{\text{*l}}}} = {\sigma _{{\text{*r}}}} = \sigma _{\text{*}}^{{\text{trial}}}$ 和速度相等 ${u_{{\text{*l}}}} = {u_{{\text{*r}}}} = u_{\text{*}}^{{\text{trial}}}$ ，按照HLLC构造格式^[14]，可以推出仅存在弹性波三波系假设下的速度和应力的HLLC黎曼解：

$\left\{\begin{gathered} \omega _{\text{l}}^{\text{E}} = - {\rho _{\text{l}}}c_{\text{l}}^{\text{E}}, \qquad\quad \omega _{\text{r}}^{\text{E}} = {\rho _{\text{r}}}c_{\text{r}}^{\text{E}} \\ u_{\text{*}}^{{\text{trial}}} = \frac{{ - \omega _{\text{l}}^{\text{E}}{u_{\text{l}}} + \omega _{\text{r}}^{\text{E}}{u_{\text{r}}} - \left( {{\sigma _{\text{l}}} - {\sigma _{\text{r}}}} \right)}}{{\omega _{\text{r}}^{\text{E}} - \omega _{\text{l}}^{\text{E}}}} \\ \sigma _{\text{*}}^{{\text{trial}}} = \frac{{ - \omega _{\text{l}}^{\text{E}}{\sigma _{\text{l}}} + \omega _{\text{r}}^{\text{E}}{\sigma _{\text{r}}} + \omega _{\text{l}}^{\text{E}}\omega _{\text{r}}^{\text{E}}\left( {{u_{\text{l}}} - {u_{\text{r}}}} \right)}}{{\omega _{\text{r}}^{\text{E}} - \omega _{\text{l}}^{\text{E}}}} \\ \rho _{{\text{*}}K}^{{\text{trial}}} = \left[ {\frac{1}{{{\rho _K}}} - \frac{{u_{\text{*}}^{{\text{trial}}} - {u_{\text{l}}}}}{{\omega _K^{\text{E}}}}} \right] \\ \end{gathered}\right.$

(7)

式中： $\omega$ 表示拉格朗日波速，下标 $K$ 表示物理量的左右状态l或者r。通过偏应力更新公式(3)和连续性方程(2)，可以得到偏应力与密度的关系：

$\frac{{{\text{d}}S}}{{{\text{d}}t}} = - \frac{{4\mu }}{3}\frac{{{\text{d}}\ln \rho }}{{{\text{d}}t}}$

(8)

积分式（8）得到接触区域的偏应力：

$S_{ {{\text{*}}{K}}}^{{\text{trial}}} = {S _{{K}}} - \frac{{4\mu }}{3}\ln {\frac{{\rho _{{{\text{*}}{K}}}^{{\text{trial}}}}}{{{\rho _{{K}}}}}}$

(9)

第二步，对所求得的偏应力进行Mises应力检测，判断接触区域有无应力屈服。如果两侧的偏应力都没达到屈服极限，那么上述应力和速度黎曼解就作为最终解代入控制方程(2)中去。如果接触波其中一侧或者两侧偏应力大于弹性极限，即接触区域产生了塑性波，则需要通过上面弹性波假设下得到的黎曼解对产生应力屈服一侧的密度、速度以及应力进行重新更新。偏应力更新公式如下：

${S_{ {{\text{*}}{K}}}} = \frac{2}{3}{Y_0}\frac{{\left| {S_{ {{\text{*}}{K}}}^{{\text{trial}}}} \right|}}{{S_{ {{\text{*}}{K}}}^{{\text{trial}}}}}$

(10)

通过再次对公式(9)进行积分，得到产生塑性变化一侧的密度：

${\rho _{{{\text{*}}{K}}}} = {\rho _{{K}}}{{\mathrm{e}}^{ - \tfrac{{4\mu }}{3}\left( {{S_{{{\text{*}}{K}}}} - {S_{{K}}}} \right)}}$

(11)

然后，通过跨越弹性波的Rankine-Hugoniot条件和密度，可以计算出速度，进而计算得到应力：

$\begin{gathered} \omega _{\text{l}}^{\text{P}} = - {\rho _{\text{l}}}c_{\text{l}}^{\text{P}}, \qquad\quad \omega _{\text{r}}^{\text{P}} = {\rho _{\text{r}}}c_{\text{r}}^{\text{P}} \\ u_{{K}}^{\text{P}} = {u_{\text{r}}} + \omega _{{K}}^{\text{P}}\left( {\frac{1}{{{\rho _{{K}}}}} - \frac{1}{{{\rho _{{{\text{*}}{K}}}}}}} \right) \\ \sigma _{{K}}^{\text{P}} = {\sigma _{{K}}} - \omega _{{K}}^{\text{P}}\left( {u_{{K}}^{\text{P}} - {u_{{K}}}} \right) \\ \end{gathered}$

(12)

在计算第三步之前，为了计算简便，先对上面已有的左右侧数据进行整合。对于接触波左侧或者右侧状态，如果该侧在第二步检测中并未达到屈服极限，则该侧的状态（上标W）为弹性状态E，即：

$\omega _K^W = \omega _K^{\text{E}}, \qquad\quad \sigma _K^W = {\sigma _K}, \qquad\quad u_K^W = {u_K}$

(13)

如果该侧到达了屈服极限，则该侧状态为塑性状态P：

$\omega _K^W = \omega _K^{\text{P}}, \qquad\quad \sigma _K^W = \sigma _K^{\text{P}}, \qquad\quad u_K^W = u_K^{\text{P}}$

(14)

第三步，通过上面得到的数据再计算一次应力和速度的黎曼解。对于新的左右状态再一次求解，得到最后的速度和应力HLLC黎曼解：

$\begin{gathered} {u_*} = \frac{{ - \omega _{\text{l}}^Wu_{\text{l}}^W + \omega _{\text{r}}^W u_{\text{r}}^W - \left( {\sigma _{\text{l}}^W - \sigma _{\text{r}}^W} \right)}}{{\omega _{\text{r}}^W - \omega _{\text{l}}^W}} \\ {\sigma _*} = \frac{{ - \omega _{\text{l}}^W\sigma _{\text{l}}^W + \omega _{\text{r}}^W\sigma _{\text{r}}^W + \omega _{\text{l}}^W\omega _{\text{r}}^W\left( {u_{\text{l}}^W - u_{\text{r}}^W} \right)}}{{\omega _{\text{r}}^W - \omega _{\text{l}}^W}} \\ \end{gathered}$

(15)

2. 算例验证

2.1 传统SPH方法

传统一维SPH方法的常用弱形式离散化控制方程如下：

$\left\{ \begin{gathered} \frac{{{\text{d}}{\rho _i}}}{{{\text{d}}t}} = {m_j}\left( {{\boldsymbol{u}}_i^\beta - {\boldsymbol{u}}_j^\beta } \right){\boldsymbol{L}}_i^{\beta \gamma }\frac{{\partial {W_{ij}}}}{{\partial {\boldsymbol{x}}_i^\gamma }} \\ \frac{{{\text{d}}{{\boldsymbol{u}}_i}}}{{{\text{d}}t}} = \sum\limits_{j = 1}^n {{m_j}\left( {\frac{{{\boldsymbol{\sigma }}_i^{\alpha \beta }}}{{{\rho _i}^2}} + \frac{{{\boldsymbol{\sigma }}_j^{\alpha \beta }}}{{{\rho _j}^2}} - {\varPi _{ij}}{{\boldsymbol{I}}^{\alpha \beta }}} \right){\boldsymbol{L}}_i^{\beta \gamma }\frac{{\partial {W_{ij}}}}{{\partial {\boldsymbol{x}}_i^\gamma }}} \\ \frac{{{\text{d}}{e_i}}}{{{\text{d}}t}} = - \frac{1}{2}\sum\limits_{j = 1}^n {{m_j}\left( {{\boldsymbol{u}}_j^\alpha - {\boldsymbol{u}}_i^\alpha } \right)\left( {\frac{{{\boldsymbol{\sigma }}_i^{\alpha \beta }}}{{{\rho _i}^2}} + \frac{{{\boldsymbol{\sigma }}_j^{\alpha \beta }}}{{{\rho _j}^2}} - {\varPi _{ij}}{{\boldsymbol{I}}^{\alpha \beta }}} \right){\boldsymbol{L}}_i^{\beta \gamma }\frac{{\partial {W_{ij}}}}{{\partial {\boldsymbol{x}}_i^\gamma }}} \\ {\boldsymbol{L}}_i^{\beta \gamma } = - {\left( {\sum\limits_{j = 1}^n {{\boldsymbol{x}}_{ij}^\beta \frac{{\partial {W_{ij}}}}{{\partial {\boldsymbol{x}}_i^\gamma }}} {V_j}} \right)^{ - 1}} \\ \end{gathered} \right.$

(16)

式中：人工黏性项 ${\varPi _{ij}}$ 在传统SPH算法中起到了稳定算法、提高计算准确性的作用，其表达方程为：

${\varPi }_{ij}=\left\{ \begin{array}{ll}\dfrac{-\alpha {\overline{C}}_{ij}{\pi}_{ij}+\beta {\pi}_{ij}^{2}}{{\overline{\rho }}_{ij}}&\qquad {u}_{ij}{x}_{ij} {\text{＜}} 0\\ 0& \qquad 其他\end{array}\right.$

(17)

式中： $\alpha$ 和 $\beta$ 为黏性力控制系数，在本文中两者取相同值，其值分别取0.5、1.0、2.0进行计算，并将结果与本文SPH-HLLC算法进行比对。另外，为了与传统计算方法一致，本文中的黎曼算法也同样在 ${u_{ij}}{x_{ij}} {\text{＜}}0$ 时采用，而当 ${u_{ij}}{x_{ij}} {\text{≥}} 0$ 时，式(2)中的黎曼解由左右状态的算术平均值 $\bar \sigma = {{\left( {{\sigma _i} + {\sigma _j}} \right)} / 2}$ 替代，相当于此时回归了无人工黏性的传统算法。

2.2 一维弹塑性数值模拟

首先，对铝的一维固体碰撞问题进行数值模拟以检测算法准确性，其初始条件如表1所示。

表 1 弹塑性碰撞测试

Table 1. Parameters related Al impact

算例		ρ/(kg·m⁻³)	u/(m·s⁻¹)	p/MPa	S/MPa	坐标区间/m	t/ms
铝碰撞	左侧	2700	200	0	−200	−1～0	0.1
铝碰撞	右侧	2700	−200	0	0	0～1	0.1

下载: 导出CSV

| 显示表格

压力通过采用Stiffened-Gas 状态方程来更新：

$p\left( {\rho ,e} \right) = c_0^2\left( {\rho - {\rho _0}} \right) + \left( {\varGamma - 1} \right)\rho e$

(18)

式中：铝的相关物理参数为初始密度 ${\rho _0} = 2\;700\;{{{\text{kg}}}/ {{{\text{m}}^{\text{3}}}}}$ ，初始声速 ${c_0} = 5\;380\;{{\text{m}}/{\text{s}}}$ ，Г=2.67，弹性模量 $\mu = 26.5\;{\text{GPa}}$ ，屈服极限 ${Y_0} = 300\;{\text{MPa}}$ 。可以看到，该算例左侧偏应力初始状态刚达到屈服极限，而右侧偏应力为零。计算结果如图1所示。

图 1 铝碰撞算例的密度、压强、速度和偏应力计算结果

Figure 1. Density, pressure, velocity and deviatoric stress profile of Al collision test

下载: 全尺寸图片幻灯片

参考解来自Chen等^[15]的2000网格计算结果。可以看出，在碰撞算例中，SPH-HLLC黎曼求解器求得的结果与传统方法黏性参数 $\alpha ,\beta = 1.0$ 接近，两种方法都能较好地能捕捉到弹塑性激波。

第二个和第三个算例分别为铝的Wilkins算例^[16]和铜铝材料撞击算例，其中参考解为采用4000网格的黎曼算法^[8,17]求得。

压强通过Mie-Grüneisen状态方程更新，其表达式为：

$p\left( {\rho ,e} \right) = {\rho _0}a_0^2f\left( \eta \right) + {\rho _0}{\varGamma _0}e$

(19)

式中： $f\left( \eta \right) = \dfrac{{\left( {\eta - 1} \right)[ {\eta - {\varGamma _0}{{( {\eta - 1} )} / 2}} ]}}{{{{[ {\eta - s( {\eta - 1} )} ]}^2}}}$ ， $\eta = \dfrac{\rho }{{{\rho _0}}}$ ， ${\rho _0}$ 、 ${a_0}$ 、 $s$ 、 ${\varGamma _0}$ 都是Mie-Grüneisen状态方程的材料常数。Wilkins算例初始状态和Mie-Grüneisen状态方程的材料常数在表2中已经列出，其中压强和内能的初始值都为0，算例的物理过程持续到时刻 $t$ 计算结束，粒子数为1000个，呈均匀分布。计算结果如图2所示。

表 2 Wilkins算例的相关参数

Table 2. Parameters related to Wilkins test

算例		ρ/(kg·m⁻³)	u/(m·s⁻¹)	Y₀/MPa	μ/GPa	a₀/(m·s⁻¹)	ρ₀/(kg·m⁻³)	s	Γ₀	坐标区间/mm	t/μs
Wilkins	左侧	2785	800	300	27.6	5328	2785	1.338	2	0 ～ 5	5
Wilkins	右侧	2785	0	300	27.6	5328	2785	1.338	2	5 ～ 50	5

下载: 导出CSV

| 显示表格

图 2 Wilkins算例的密度、压强、速度和偏应力计算结果

Figure 2. Density, pressure, velocity and deviatoric stress profile of Wilkins test

下载: 全尺寸图片幻灯片

从图2中可以看出，SPH-HLLC黎曼求解器耦合算法能较准确地进行一维模拟。在 $x=0.032$ 左侧这一段是反射稀疏波作用区域，而在右侧 $0.032\text{＜}x\text{＜}0.04$ 区域可以看到右行的弹性波和塑性波。可以看到，传统算法中，当人工黏性参数 $\alpha ,\beta = 0.5$ 时，塑性激波波后出现了震荡，且对于塑性稀疏波的捕捉出现了失真。对于右侧激波，SPH-HLLC算法的准确性与传统方法取人工黏性参数 $\alpha ,\beta \in [1.0,2.0]$ 的结果接近。

第二个算例是铜铝材料固体碰撞算例，初始状态与Mie-Grüneisen状态方程的材料常数在表3中已经列出，其中压强和内能的初始值都为零，粒子数为1000个，呈均匀分布。计算结果如图3所示。

表 3 铜铝材料撞击算例的相关参数

Table 3. Parameters related to Cu/Al collision

算例		ρ/(kg·m⁻³)	u/(m·s⁻¹)	Y₀/MPa	μ/GPa	a₀/(m·s⁻¹)	ρ₀/(kg·m⁻³)	s	Γ₀	坐标区间/mm	t/μs
铜铝材料撞击	左侧	8930	60	90	45.0	3940	8930	1.490	2	0 ～ 25	2
铜铝材料撞击	右侧	2785	0	300	27.6	5328	2785	1.338	2	25～50	2

下载: 导出CSV

| 显示表格

图 3 铜铝撞击算例的压强、速度和偏应力计算结果

Figure 3. Pressure, velocity and deviatoric stress profile of Cu/Al collision test

下载: 全尺寸图片幻灯片

可以看到，在激波捕捉精度上，SPH-HLLC算法的准确性与人工黏性参数取 $\alpha ,\beta \in [1.0,2.0]$ 时得到的结果较为接近。而对于两种不同材料的接触界面，常用的传统弱形式算法很明显出现了巨大的压强和偏应力的震荡，且将人工黏性参数从0.5增大到2.0对减少震荡没有产生显著效果，而新算法显著地抑制了接触界面的非物理震荡。其物理机制在于，新算法通过构造黎曼解格式沿特征波方向加入了隐式黏性，使得压强和偏应力的震荡得到了较好的抑制。而传统人工黏性则只能沿着粒子速度靠近方向增加显式黏性，无法区分特征波方向，强度依靠统一经验参数设定，只能靠趋近速度微调，因此，对复杂波系的非物理震荡抑制效果有限。在计算效率方面，传统方法以 $\alpha ,\beta = 1$ 为例与SPH-HLLC算法对比，Wilkins算例和铜铝材料碰撞算例的计算时长分别为13.45、15.13 s和23.71、27.16 s。可以看到，SPH-HLLC算法相较于传统方法计算时间略有增加。

2.3 二维厚壁圆筒铍壳体碰撞测试

该算例最初由Howell等^[18]提出。它计算了一个全场拥有指向圆心的轴向初速度的厚壁圆筒铍壳体的动力学过程，其中初动能会在计算结束时通过塑性耗散完全转化为内能。其状态方程使用Mie-Grüneisen模型，物理参数如表4所示。

表 4 二维厚壁圆筒铍壳体碰撞算例的相关参数

Table 4. Parameters related to collapse of a thick-walled cylindrical beryllium shell

算例	ρ₀/(kg·m⁻³)	\|u\|/(m·s⁻¹)	Y₀/MPa	μ/GPa	a₀/(m·s⁻¹)	s	Γ	t/μs
铍壳体	1845	417.1	330	151.9	12870	1.124	2	130

下载: 导出CSV

| 显示表格

对粒子初始直径 $l_0=0.2,\ 0.125,\ 0.08\; \text{mm}$ 的三种粒度分别用黏性参数 $\alpha,\ \beta=1.0$ 的传统方法以及新SPH-HLLC算法计算，共计 $3 \times 2 = 6$ 个算例。

由于该算例是二维算例，需要对每个粒子对单独进行上述的一维HLLC黎曼求解。首先，要得到每个粒子对ij的黎曼间断初始条件：

${\sigma _K} = {{\boldsymbol{n}}_{ij}} \cdot {{\boldsymbol{\sigma }}_K} \cdot {{\boldsymbol{n}}_{ij}},\qquad\quad {{\boldsymbol{S}}_K} = {{\boldsymbol{n}}_{ij}} \cdot {{\boldsymbol{S}}_K} \cdot {{\boldsymbol{n}}_{ij}},\qquad\quad {u_K} = {{\boldsymbol{u}}_K} \cdot {{\boldsymbol{n}}_{ij}}$

(20)

式中： ${\boldsymbol{\sigma }}$ 和 ${\boldsymbol{S}}$ 为应力和偏应力张量， ${{\boldsymbol{n}}_{ij}} = \dfrac{{{x_{ij}}}}{{| {{x_{ij}}} |}}$ 为由i指向j的单位向量。将上述粒子对ij的一维黎曼间断初始条件代入2.2节中进行计算，就可以求得一维方向的应力黎曼解 ${\sigma _*}$ ，然后通过下式得到应力黎曼解张量 ${{\boldsymbol{\sigma }}_*}$ ：

$\begin{gathered}\boldsymbol{\sigma}_{\mathrm{\mathit{K},pa}}=\boldsymbol{\sigma}_K-\sigma_K\left(\boldsymbol{n}_{ij}\otimes\boldsymbol{n}_{ij}\right) \\ \boldsymbol{\sigma}_*=\frac{1}{2}\left(\boldsymbol{\sigma}_{\mathrm{l,pa}}+\boldsymbol{\sigma}_{\mathrm{r,pa}}\right)+\sigma_*\left(\boldsymbol{n}_{ij}\otimes\boldsymbol{n}_{ij}\right) \\ \end{gathered}$

(21)

该算例的初始分布以及计算结果空间分布示意图如图4所示。

图 4 二维厚壁圆筒铍壳体碰撞测试的初始分布、传统算法结果以及新SPH-HLLC算法结果

Figure 4. Initial distribution, traditional algorithm results, and new SPH-HLLC algorithm results of two-dimensional thick walled cylindrical beryllium shell collision test

下载: 全尺寸图片幻灯片

可以看到，传统方法和新SPH-HLLC算法计算得到的结果和理论解^[19]都与计算的理论内径0.05 m和外径0.0781 m较为接近，说明两种算法在二维情况下都有较强的适用性。当 ${l_0} = 0.08\;{\text{mm}}$ 时，两种方法的计算时长分别为941、1242 s，SPH-HLLC方法计算时长是传统方法计算时长的约1.3倍。当然，事实上，计算时长受到很多因素的影响，区域划分、并行效率等都会影响到采用不同方法的计算时长，从这里可以粗略地看出本文方法并不会显著增加计算负担。

图5为不同分辨率下传统方法和新SPH-HLLC算法的无量纲能量E^d（能量除以初始总能）随无量纲时间t^d（时间除以模拟停止时间）的演化图。图中 $E_k^{\mathrm{d}}$ 为无量纲动能， $E_i^{\mathrm{d}}$ 为无量纲内能， $E_t^{\mathrm{d}}$ 为无量纲总能。

图 5 二维厚壁圆筒铍壳体碰撞测试的无量纲能量随无量纲时间的变化

Figure 5. Variation of non-dimensional energy with non-dimensional time in collision test of two-dimensional thick walled cylindrical beryllium shells

下载: 全尺寸图片幻灯片

从图5中可以看到，在动能转化上，传统方法和新SPH-HLLC算法的动能随时间变化曲线随着粒子分辨率的提高是渐进接近精确解的。但传统方法的内能略低于精确解，且随粒子间距的增大而增大，而新SPH-HLLC算法的结果是内能略高于精确解，在求解最后却是一致趋近于精确解，与动能变化一致。这说明新SPH-HLLC算法在求解过程中由于自适应调节了数值黏性，动能很好地转化为了内能，而不是如传统算法将耗散转换成熵增，导致部分内能损失。另外，注意到传统算法和新SPH-HLLC算法都出现了一定程度的总能不守恒问题，这是因为，我们的控制方程使用修复矩阵 ${\boldsymbol{L}}$ ，该算法修复了自由界面的截断误差^[20-21]，但破坏了SPH构造格式的对称性，因此带来了一定的能量不守恒问题，表5定量给出了两种方法结束时刻与初始时刻的总能对比，可见，新SPH-HLLC算法较好地控制了守恒偏差不超过1%，而传统算法由于熵增耗散，导致守恒偏差超过5%，即使更细密的分辨率也难以弥补。

表 5 计算前后总能比（结束时刻能量/初始时刻能量）

Table 5. Total energy ratio (final total energy/initial total energy)

方法	总能比/%
方法	l₀ = 0.2 mm	l₀ = 0.125 mm	l₀ = 0.08 mm
传统	90.40	92.75	94.69
SPH-HLLC	100.59	100.72	100.64

下载: 导出CSV

| 显示表格

3. 结论与展望

提出了一种一维理想弹塑性条件下构造简单的SPH-HLLC方法，并用使用了不同状态方程的三个一维碰撞算例和一个二维算例对新方法的准确性进行了验证，且与传统SPH方法进行了对比。计算结果表明：

(1)本文构造的SPH-HLLC算法能较准确地对一维完美弹塑性材料碰撞进行模拟，其对于激波捕捉的精度接近传统SPH方法（人工黏性项参数 $\alpha ,\beta \in [1.0,2.0]$ ）；

(2)本文构造的SPH-HLLC算法在不同材料接触界面展现出了比传统SPH方法更优异的压强及偏应力非物理震荡抑制效果，且新的SPH-HLLC算法相比传统方法在总能守恒上具有更好的表现；

(3)通过在二维计算中对每一个粒子对进行一维黎曼求解，本文将该一维黎曼解推广到二维算例中，并且取得了正确、合理的结果。

需要指出的是，本文提出的方法可以进一步推广至一般弹塑性模型，但需要处理后继屈服极限随应力路径变化的问题

图 1 铝碰撞算例的密度、压强、速度和偏应力计算结果

Figure 1. Density, pressure, velocity and deviatoric stress profile of Al collision test

下载: 全尺寸图片幻灯片

图 2 Wilkins算例的密度、压强、速度和偏应力计算结果

Figure 2. Density, pressure, velocity and deviatoric stress profile of Wilkins test

下载: 全尺寸图片幻灯片

图 3 铜铝撞击算例的压强、速度和偏应力计算结果

Figure 3. Pressure, velocity and deviatoric stress profile of Cu/Al collision test

下载: 全尺寸图片幻灯片

图 4 二维厚壁圆筒铍壳体碰撞测试的初始分布、传统算法结果以及新SPH-HLLC算法结果

Figure 4. Initial distribution, traditional algorithm results, and new SPH-HLLC algorithm results of two-dimensional thick walled cylindrical beryllium shell collision test

下载: 全尺寸图片幻灯片

图 5 二维厚壁圆筒铍壳体碰撞测试的无量纲能量随无量纲时间的变化

Figure 5. Variation of non-dimensional energy with non-dimensional time in collision test of two-dimensional thick walled cylindrical beryllium shells

下载: 全尺寸图片幻灯片

表 1 弹塑性碰撞测试

Table 1. Parameters related Al impact

算例		ρ/(kg·m⁻³)	u/(m·s⁻¹)	p/MPa	S/MPa	坐标区间/m	t/ms
铝碰撞	左侧	2700	200	0	−200	−1～0	0.1
铝碰撞	右侧	2700	−200	0	0	0～1	0.1

下载: 导出CSV

表 2 Wilkins算例的相关参数

Table 2. Parameters related to Wilkins test

算例		ρ/(kg·m⁻³)	u/(m·s⁻¹)	Y₀/MPa	μ/GPa	a₀/(m·s⁻¹)	ρ₀/(kg·m⁻³)	s	Γ₀	坐标区间/mm	t/μs
Wilkins	左侧	2785	800	300	27.6	5328	2785	1.338	2	0 ～ 5	5
Wilkins	右侧	2785	0	300	27.6	5328	2785	1.338	2	5 ～ 50	5

下载: 导出CSV

表 3 铜铝材料撞击算例的相关参数

Table 3. Parameters related to Cu/Al collision

算例		ρ/(kg·m⁻³)	u/(m·s⁻¹)	Y₀/MPa	μ/GPa	a₀/(m·s⁻¹)	ρ₀/(kg·m⁻³)	s	Γ₀	坐标区间/mm	t/μs
铜铝材料撞击	左侧	8930	60	90	45.0	3940	8930	1.490	2	0 ～ 25	2
铜铝材料撞击	右侧	2785	0	300	27.6	5328	2785	1.338	2	25～50	2

下载: 导出CSV

表 4 二维厚壁圆筒铍壳体碰撞算例的相关参数

Table 4. Parameters related to collapse of a thick-walled cylindrical beryllium shell

算例	ρ₀/(kg·m⁻³)	\|u\|/(m·s⁻¹)	Y₀/MPa	μ/GPa	a₀/(m·s⁻¹)	s	Γ	t/μs
铍壳体	1845	417.1	330	151.9	12870	1.124	2	130

下载: 导出CSV

表 5 计算前后总能比（结束时刻能量/初始时刻能量）

Table 5. Total energy ratio (final total energy/initial total energy)

方法	总能比/%
方法	l₀ = 0.2 mm	l₀ = 0.125 mm	l₀ = 0.08 mm
传统	90.40	92.75	94.69
SPH-HLLC	100.59	100.72	100.64

下载: 导出CSV

参考文献(21)

[1]	MONAGHAN J J. Smoothed particle hydrodynamics [J]. Annual Review of Astronomy and Astrophysics, 1992, 30: 543–574. DOI: 10.1146/annurev.aa.30.090192.002551.
[2]	VILA J P. On particle weighted methods and smooth particle hydrodynamics [J]. Mathematical Models and Methods in Applied Sciences, 1999, 9(2): 161–209. DOI: 10.1142/S0218202599000117.
[3]	PARSHIKOV A N, MEDIN S A. Smoothed particle hydrodynamics usinginterparticle contact algorithms [J]. Journal of Computational Physics, 2002, 180(1): 358–382. DOI: 10.1006/jcph.2002.7099.
[4]	LIBERSKY L D, RANDLES P W. Shocks and discontinuities in particle methods [J]. AIP Conference Proceedings, 2006, 845(1): 1089–1092. DOI: 10.1063/1.2263512.
[5]	MEHRA V, CHATURVEDI S. High velocity impact of metal sphere on thin metallic plates: a comparative smooth particle hydrodynamics study [J]. Journal of Computational Physics, 2006, 212(1): 318–337. DOI: 10.1016/j.jcp.2005.06.020.
[6]	LIN X, BALLMANN J. A Riemann solver and a second-order Godunov method for elastic-plastic wave propagation in solids [J]. International Journal of Impact Engineering, 1993, 13(3): 463–478. DOI: 10.1016/0734-743X(93)90118-Q.
[7]	姚成宝, 付梅艳, 韩峰, 等. 基于多介质Riemann问题的流体-固体耦合数值方法及其在爆炸与冲击问题中的应用 [J]. 兵工学报, 2021, 42(2): 340–355. DOI: 10.3969/j.issn.1000-1093.2021.02.012. YAO C B, FU M Y, HAN F, et al. A numerical scheme for fluid-solid interactions based on multi-medium Riemann problem and its application in explosion andimpact problems [J]. Acta Armamentarii, 2021, 42(2): 340–355. DOI: 10.3969/j.issn.1000-1093.2021.02.012.
[8]	CHENG J B. Harten-Lax-van Leer-contact (HLLC) approximation Riemann solver with elastic waves for one-dimensional elastic-plastic problems [J]. Applied Mathematics and Mechanics, 2016, 37(11): 1517–1538. DOI: 10.1007/s10483-016-2104-9.
[9]	GAO S, LIU T G. 1D exact elastic-perfectly plastic solid Riemann solver and its multi-material application [J]. Advances in Applied Mathematics and Mechanics, 2017, 9(3): 621–650. DOI: 10.4208/aamm.2015.m1340.
[10]	GAO S, LIU T G, YAO C B. A complete list of exact solutions for one-dimensional elastic-perfectly plastic solid Riemann problem without vacuum [J]. Communications in Nonlinear Science and Numerical Simulation, 2018, 63: 205–227. DOI: 10.1016/j.cnsns.2018.02.030.
[11]	LI X, ZHAI J Y, SHEN Z J. An HLLC-type approximate Riemann solver for two-dimensional elastic-perfectly plastic model [J]. Journal of Computational Physics, 2022, 448: 110675. DOI: 10.1016/j.jcp.2021.110675.
[12]	LIU M B, LIU G R. Smoothed particle hydrodynamics (SPH): an overview and recent developments [J]. Archives of Computational Methods in Engineering, 2010, 17(1): 25–76. DOI: 10.1007/s11831-010-9040-7.
[13]	HUI W H, KUDRIAKOV S. On wall overheating and other computational difficulties of shock-capturing methods [J]. Computational Fluid Dynamics Journal, 2001, 10(2): 192–209.
[14]	TORO E F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction [M]. New York: Springer, 1997.
[15]	CHEN Q, LI L, QI J, et al. A cell-centeredLagrangian scheme with an elastic-perfectly plastic solid Riemann solver for wave propagations in solids [J]. Advances in Applied Mathematics and Mechanics, 2022, 14(3): 703–724. DOI: 10.4208/aamm.OA-2020-0344.
[16]	WILKINS M L. Calculation of elastic-plastic flow: NSA-18-002406 [R]. Livermore: Lawrence Radiation Laboratory, 1963.
[17]	LIU L, CHENG J B, LIU Z. A multi-material HLLC Riemann solver with both elastic and plastic waves for 1D elastic-plastic flows [J]. Computers & Fluids, 2019, 192: 104265. DOI: 10.1016/j.compfluid.2019.104265.
[18]	HOWELL B P, BALL G J. A free-Lagrange augmented Godunov method for the simulation of elastic-plastic solids [J]. Journal of Computational Physics, 2002, 175(1): 128–167. DOI: 10.1006/jcph.2001.6931.
[19]	MAIRE P H, ABGRALL R, BREIL J, et al. A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two-dimensional unstructured grids [J]. Journal of Computational Physics, 2013, 235: 626–665. DOI: 10.1016/j.jcp.2012.10.017.
[20]	QUINLAN N J, BASA M, LASTIWKA M. Truncation error in mesh-free particle methods [J]. International Journal for Numerical Methods in Engineering, 2006, 66(13): 2064–2085. DOI: 10.1002/nme.1617.
[21]	ZHANG Z L, LIU M B. Smoothed particle hydrodynamics with kernel gradient correction for modeling high velocity impact in two- and three-dimensional spaces [J]. Engineering Analysis with Boundary Elements, 2017, 83: 141–157. DOI: 10.1016/j.enganabound.2017.07.015.

施引文献

资源附件(0)

访问统计

图(5) / 表(5)

计量

文章访问数: 302
HTML全文浏览量: 93
PDF下载量: 90
被引次数: 0

1. 计算方法
1.1 SPH控制方程及其离散化
1.2 一维SPH-HLLC黎曼解
2. 算例验证
2.1 传统SPH方法
2.2 一维弹塑性数值模拟
2.3 二维厚壁圆筒铍壳体碰撞测试
3. 结论与展望

一维理想弹塑性体的SPH-HLLC耦合算法

doi: 10.11883/bzycj-2024-0004

作者简介: 王展铭（1997－ ），男，博士研究生，wangzmrr@mail.ustc.edu.cn

通讯作者: 黄生洪（1974－ ），男，博士，研究员，博士生导师，hshnpu@ustc.edu.cn

计量

出版历程

SPH-HLLC coupled method for one-dimentional elastic-perfectly plastic model

1. 计算方法

1.1 SPH控制方程及其离散化

1.2 一维SPH-HLLC黎曼解

2. 算例验证

2.1 传统SPH方法

2.2 一维弹塑性数值模拟

2.3 二维厚壁圆筒铍壳体碰撞测试

3. 结论与展望

计量

出版历程

目录

1. 计算方法

1.1 SPH控制方程及其离散化

1.2 一维SPH-HLLC黎曼解

2. 算例验证

2.1 传统SPH方法

2.2 一维弹塑性数值模拟

2.3 二维厚壁圆筒铍壳体碰撞测试

3. 结论与展望

作者简介:
王展铭（1997－　），男，博士研究生，wangzmrr@mail.ustc.edu.cn

通讯作者:
黄生洪（1974－　），男，博士，研究员，博士生导师，hshnpu@ustc.edu.cn