一种在欧拉坐标系下计算一维理想弹塑性固体的技术

摘要

本发明提出了一种在欧拉坐标系下计算一维理想弹塑性固体的技术，其发明的主要内容为一套完整的在欧拉坐标系下利用一维理想弹塑性固体的模型计算其相关物理量的技术。本发明的创新点主要体现在一维Hooke定律中的物质导数在欧拉坐标系下的计算方式。本发明的提出，可直接用来计算一维理想弹塑性固体，并在一维理想弹塑性固体受外力作用、与其他介质耦合作用等实际工程应用中具有重要意义。

说明书

技术领域

本发明涉及一种计算一维理想弹塑性固体的技术，具体涉及一种在欧拉坐标系下计算一维理想弹塑性 固体的技术。

背景技术

一维理想弹塑性固体模型可以较准确描述固体(如铝、钢等金属)在受一般强度的外力作用下的各物 理量变化情况。因此，研究理想弹塑性固体的计算技术，在实际工程中具有重要的应用价值和广泛的应用 前景。

目前，虽然已经存在一些计算理想弹塑性固体的技术，但均与本发明提出的技术不同。比如， M.L.Wilkins在1964年在提出理想弹塑性固体的模型之后，采用有限差分法对该模型进行求解，其中使用 了复杂的全离散形式。再如，B.P.Howell在2000年采用Free Lagrange方法对理想弹塑性固体进行计 算。该方法是在拉格朗日坐标下进行计算的，虽然在计算一些变量(如偏应力)上可以得到简化，但推广 到高维时变得十分复杂。为了使理想弹塑性固体的计算既简单有准确，本发明直接在欧拉坐标系下进行计 算，只需将Hooke定律中的导数作为物质导数进行处理。值得一提的是，本发明是在1993年M.B.Tyndall 的工作上受到启发而提出的。然而，M.B.Tyndall的计算方法中却存在一些错误。首先，在计算欧拉坐 标下每个网格点在上一时间步的位置时，他在固定的网格点上对速度采取了时间平均。这种计算方法在拉 格朗日坐标下是正确的，在欧拉坐标下却是错误的，与他建立的欧拉坐标下的控制方程相矛盾。其次，他 在计算相关物理量在每个网格点上一时间步处的值的时候，采用了该点两侧网格点的抛物插值(二次函数 插值)。可是，理想弹塑性固体的控制方程是双曲方程，传播具有方向性，采用抛物插值会造成计算的不 准确，甚至造成不稳定并产生错误。实际的数值计算也验证了他的方法确实存在一些错误。对于这一问题， 本发明则直接采取迎风线性插值。总之，本发明提出的一维理想弹塑性固体的计算技术兼顾了方法的简单 性与正确性。

发明内容

本发明提出的的计算一维理想弹塑性固体的技术，其发明内容主要体现在在欧拉坐标系下的一套完整 的计算一维理想弹塑性固体的技术，其创新点主要体现在一维Hooke定律中的物质导数在欧拉坐标系下的 计算方式。

对于一维情况，理想弹塑性固体在欧拉坐标系下的控制方程为

$\frac{\partial U}{\partial t}+\frac{\partial F\left(U\right)}{\partial x}=0---\left(1\right)$

$U=\left(\begin{array}{c}\rho \\ \mathrm{\rho u}\\ E\end{array}\right),F\left(U\right)=\left(\begin{array}{c}\mathrm{\rho u}\\ {\mathrm{\rho u}}^2-{\sigma }_x\\ (E-{\sigma }_x)u\end{array}\right)$

此处，ρ是密度，u是速度，p是压力，E是总能，σ_x是x方向的总应力。此外，对于理想弹塑性固体， 其总应力和压力还满足下面关系：

σ_x＝-p+s_x

其中，s_x是x方向的偏应力。当理想弹塑性固体处于弹性状态，有

$\stackrel{·}{p}=K\frac{\stackrel{·}{\rho }}{\rho }$和${\stackrel{·}{s}}_x=-\frac{4}{3}\mu \frac{\stackrel{·}{\rho }}{\rho }$

其中K是体积模量，μ是剪切模量。当理想弹塑性固体处于塑性状态，有

$p=c_0^2(\rho -{\rho }_0)+({\gamma }_s-1)\mathrm{\rho e}$和$s_x=\pm \frac{2}{3}{Y}_0$

其中c₀，ρ₀，γ_s均为与具体固体有关的常数，Y₀是屈服强度；对于偏应力s_x，正号表示固体处于拉伸状 态，符号表示固体处于压缩状态。当理想弹塑性固体满足下面方程时为弹性状态

${s_x}^2,,{\left(\frac{2}{3}{Y}_0\right)}^2$

当上述不等式不成立时，固体处于塑性状态。

本发明的具体发明内容可以归结为如下计算技术。假设已知一维理想弹塑性固体在第n个时间步的各 变量值需要将这些变量值推进到第n+1个时间步，得到其计算技术通过以下六个步骤来实现：

1.求解控制方程(1)，将控制方程中第n个时间步的各变量值更新至第n+1个 时间步，得到${\rho }_i^{n+1},u_i^{n+1},E_i^{n+1}.$

2.计算欧拉坐标下每个网格点在第n个时间步的位置，记为x_old，有

3.采用迎风线性插值，计算ρ，p，s_x在x_old处的值，记作ρ_old，p_old，s_xold，如

$\left(\begin{array}{c}{\rho }_{\mathrm{old}}={\rho }_i^n-\frac{{\rho }_i^n-{\rho }_{i-1}^n}{\mathrm{\Delta x}}(x_i^{n+1}-x_{\mathrm{old}}),u_i^{n+1}\ge 0\\ {\rho }_{\mathrm{old}}={\rho }_i^n+\frac{{\rho }_{i+1}^n-{\rho }_i^n}{\mathrm{\Delta x}}(x_{\mathrm{old}}-x_i^{n+1}),u_i^{n+1}<0\end{array}\right)$

p_old和s_xold也可通过类似方式进行计算。

4.利用Hooke定律及迎风线性插值，得到初步的

$\left(\begin{array}{c}{s}_{\mathrm{xi}}^{n+1}=s_{\mathrm{xold}}+2\mu [\frac{u_i^{n+1}-u_{i-1}^{n+1}}{\mathrm{\Delta x}}\mathrm{\Delta t}+\frac{1}{3}\ln \left(\frac{{\rho }_i^{n+1}}{{\rho }_{\mathrm{old}}}\right)],u_i^{n+1}\ge 0\\ s_{\mathrm{xi}}^{n+1}=s_{\mathrm{xold}}+2\mu [\frac{u_{i+1}^{n+1}-u_i^{n+1}}{\mathrm{\Delta x}}\mathrm{\Delta t}+\frac{1}{3}\ln \left(\frac{{\rho }_i^{n+1}}{{\rho }_{\mathrm{old}}}\right)],u_i^{n+1}<0\end{array}\right)$

5.通过von Mises屈服条件判断在每一欧拉网格点处的弹塑性状态，并更新压力值至如 果某一网格点处满足von Mises屈服条件，即则固体处于弹性状态，并且压力通过 Hooke定律计算

$p_i^{n+1}=p_{\mathrm{old}}+K\ln \left(\frac{{\rho }_i^{n+1}}{{\rho }_{\mathrm{old}}}\right)$

若在该网格点不满足von Mises屈服条件，即则固体处于塑形状态，并且压力通过状 态方程计算

$p_i^{n+1}=c_0^2({\rho }_i^{n+1}-{\rho }_0)+({\gamma }_s-1)(E_i^{n+1}-\frac{1}{2}{\rho }_i^{n+1}{u}_i^{n+1}{u}_i^{n+1})$

同时，使偏应力满足理想塑性条件

$s_{\mathrm{xi}}^{n+1}=\frac{2}{3}{Y}_0\frac{s_{\mathrm{xi}}^{n+1}}{\left|s_{\mathrm{xi}}^{n+1}\right|}$

6.返回步骤1直到达到设定的时间迭代要求。

附图说明

图1是本发明在欧拉坐标系下计算一维理想弹塑性固体的流程图；

图2至图4是本发明在欧拉坐标系下计算一维理想弹塑性固体的一个算例结果。

具体实施方式

为了说明本发明的具体实施方式，下面将对一个算例进行演示。考虑在铝中的一维黎曼问题，其中该 黎曼问题左侧的无量纲初值为u_L＝20.0，p_L＝1.0,ρ_L＝2.7，s_L＝0.0,右侧的无量纲初值为 u_R＝-20.0,p_R＝1.0,ρ_R＝2.7,s_R＝0.0。无量纲的求解区间[0,1]中等距分布着2000个欧拉网格点， 并且黎曼问题的初始界面为0.0。同时，铝的理想弹塑性模型的相关无量纲参数分别为ρ₀＝2.71， c₀＝538.0，γ_s＝2.71，K＝740000.0，μ＝265000.0,Y₀＝3000.0。

该问题将在界面左右两侧同时产生弹性波和塑性波。取时间步长Δt＝0.0000015，采用 Lax-Friedrich格式计算，得到在时间t＝0.001时铝中的负的总应力、速度、密度如图2至图4所示。