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

Abstract

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

Description

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

技术领域

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

背景技术

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

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

同时，与一维理想弹塑性固体的计算不同的是，二维情况下模型方程均发生变化，导致其计算更加困难、计算过程更加复杂。这些都需要重新进行考虑。

发明内容

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

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

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

$U=\left(\begin{array}{c}\mathrm{\ρ}\\ \mathrm{\ρu}\\ \mathrm{\ρv}\\ E\end{array}\right),F\left(U\right)=\left(\begin{array}{c}\mathrm{\ρu}\\ \mathrm{\ρ}{u}^2-{\mathrm{\σ}}_x\\ \mathrm{\ρuv}-s_{\mathrm{xy}}\\ (E-{\mathrm{\σ}}_x)u-s_{\mathrm{xy}}v\end{array}\right),G\left(U\right)=\left(\begin{array}{c}\mathrm{\ρv}\\ \mathrm{\ρuv}-s_{\mathrm{xy}}\\ \mathrm{\ρ}{v}^2-{\mathrm{\σ}}_y\\ (E-{\mathrm{\σ}}_y)v-s_{\mathrm{xy}}u\end{array}\right)$

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

σ_x＝-p+s_x，σ_y＝-p+s_y

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

$\stackrel{\·}{\mathrm{\ρ}}=K\frac{\stackrel{\·}{\mathrm{\ρ}}}{\mathrm{\ρ}}$

${\stackrel{\·}{s}}_x=2\mathrm{\μ}({\stackrel{\·}{s}}_x+\frac{1}{3}\frac{\stackrel{\·}{\mathrm{\ρ}}}{\mathrm{\ρ}}),{\stackrel{\·}{s}}_y=\mathrm{2\μ}({\stackrel{\·}{\mathrm{\ϵ}}}_y+\frac{1}{3}\frac{\stackrel{\·}{\mathrm{\ρ}}}{\mathrm{\ρ}}),{\stackrel{\·}{s}}_{\mathrm{xy}}=\mathrm{\μ}{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{xy}}$

其中K是体积模量，μ是剪切模量，ε_x、ε_y、εx_y分别是相对应的应变，并且有

${\stackrel{\·}{\mathrm{\ϵ}}}_x=\frac{\∂u}{\∂x},{\stackrel{\·}{\mathrm{\ϵ}}}_y=\frac{\∂v}{\∂y},{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{xy}}=\frac{\∂v}{\∂x}+\frac{\∂u}{\∂y}$

当理想弹塑性固体处于塑性状态，有

$p=c_0^2(\mathrm{\ρ}-{\mathrm{\ρ}}_0)+({\mathrm{\γ}}_s-1)\mathrm{\ρe}$

$s_x=s_x\frac{\sqrt{\frac{2}{3}}{Y}_0}{\sqrt{{s_x}^2+{s_y}^2+{s_{\mathrm{xy}}}^2}},s_y=s_y\frac{\sqrt{\frac{2}{3}}{Y}_0}{\sqrt{{s_x}^2+{s_y}^2+{s_{\mathrm{xy}}}^2}},s_{\mathrm{xy}}=s_{\mathrm{xy}}\frac{\sqrt{\frac{2}{3}}{Y}_0}{\sqrt{{s_x}^2+{s_y}^2+{s_{\mathrm{xy}}}^2}}$

其中c₀，ρ₀，γ_s均为与具体固体有关的常数，Y₀是屈服强度。当理想弹塑性固体满足下面方程时为弹性状态

${s_x}^2+{s_y}^2+2{s_{\mathrm{xy}}}^2\≤{\left(\frac{2}{3}{Y}_0\right)}^2$

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

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

1.求解控制方程(1)，将控制方程中第n个时间步的各变量值更新至第n+1个时间步，得到

2.计算欧拉坐标下每个网格点()在第n个时间步的位置，记为(x_old,y_old)，有 $x_{\mathrm{old}}=x_i^{n+1}-u_{i,j}^{n+1}\mathrm{\Δt},y_{\mathrm{old}}=y_i^{n+1}-v_{i,j}^{n+1}\mathrm{\Δt}.$

3.采用迎风线性插值，计算ρ，p，s_x，s_y，s_xy在(x_old,y_old)处的值，记作ρ_old，p_old，s_xold，s_yold，s_xyold如

$\left\{\begin{array}{c}{\mathrm{\&rho;}}_{\mathrm{old}}=\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)-\frac{\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)-\mathrm{\&rho;}(x_{\mathrm{old}},y_{j-1}^n)}{\mathrm{\&Delta;y}}(y_j^n-y_{\mathrm{old}}),v_{i,j}^{n+1}\&GreaterEqual;0\\ {\mathrm{\&rho;}}_{\mathrm{old}}=\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)+\frac{\mathrm{\&rho;}(x_{\mathrm{old}},y_{j+1}^n)-\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)}{\mathrm{\&Delta;y}}(y_{\mathrm{old}}-y_j^n),v_{i,j}^{n+1}<0\end{array}\right.$

其中

$\left\{\begin{array}{c}\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)={\mathrm{\&rho;}}_{i,j}^n-\frac{{\mathrm{\&rho;}}_{i,j}^n-{\mathrm{\&rho;}}_{i-1,j}^n}{\mathrm{\&Delta;x}}(x_i-x_{\mathrm{old}})\\ \mathrm{\&rho;}(x_{\mathrm{old}},y_{j-1}^n)={\mathrm{\&rho;}}_{i,j-1}^n-\frac{{\mathrm{\&rho;}}_{i,j-1}^n-{\mathrm{\&rho;}}_{i-1,j-1}^n}{\mathrm{\&Delta;x}}(x_i-x_{\mathrm{old}}),u_{i,j}^{n+1}\&GreaterEqual;0\\ \mathrm{\&rho;}(x_{\mathrm{old}},y_{j+1}^n)={\mathrm{\&rho;}}_{i,j+1}^n-\frac{{\mathrm{\&rho;}}_{i,j+1}^n-{\mathrm{\&rho;}}_{i-1,j+1}^n}{\mathrm{\&Delta;x}}(x_i-x_{\mathrm{old}})\end{array}\right.$

$\left\{\begin{array}{c}\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)={\mathrm{\&rho;}}_{i,j}^n+\frac{{\mathrm{\&rho;}}_{i+1,j}^n-{\mathrm{\&rho;}}_{i-1,j}^n}{\mathrm{\&Delta;x}}(x_{\mathrm{old}}-x_i)\\ \mathrm{\&rho;}(x_{\mathrm{old}},y_{j-1}^n)={\mathrm{\&rho;}}_{i,j-1}^n+\frac{{\mathrm{\&rho;}}_{i+1,j-1}^n-{\mathrm{\&rho;}}_{i,j-1}^n}{\mathrm{\&Delta;x}}(x_{\mathrm{old}}-x_i),u_{i,j}^{n+1}<0\\ \mathrm{\&rho;}(x_{\mathrm{old}},y_{j+1}^n)={\mathrm{\&rho;}}_{i,j+1}^n+\frac{{\mathrm{\&rho;}}_{i+1,j+1}^n-{\mathrm{\&rho;}}_{i,j+1}^n}{\mathrm{\&Delta;x}}(x_{\mathrm{old}}-x_i)\end{array}\right.$

p_old，s_xold，s_yold，s_xyold也可通过类似方式进行计算。

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

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

类似地，也可得到

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

$p_{i,j}^{n+1}=p_{\mathrm{old}}+K\ln \left(\frac{{\mathrm{\&rho;}}_{i,j}^{n+1}}{{\mathrm{\&rho;}}_{\mathrm{old}}}\right)$

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

$p_{i,j}^{n+1}=c_0^2({\mathrm{\&rho;}}_{i,j}^{n+1}-{\mathrm{\&rho;}}_0)+({\mathrm{\&gamma;}}_s-1)(E_{i,j}^{n+1}-\frac{1}{2}{\mathrm{\&rho;}}_i^{n+1}{u}_{i,j}^{n+1}{u}_{i,j}^{n+1})$

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

$s_{\mathrm{xi},j}^{n+1}=s_{\mathrm{xi},j}^{n+1}\frac{\sqrt{\frac{2}{3}}{Y}_0}{\sqrt{s_{\mathrm{xi},j}^2+s_{\mathrm{yi},j}^2+s_{\mathrm{xyi},j}^2}},s_{\mathrm{yi},j}^{n+1}=s_{\mathrm{yi},j}^{n+1}\frac{\sqrt{\frac{2}{3}}{Y}_0}{\sqrt{s_{\mathrm{xi},j}^2+s_{\mathrm{yi},j}^2+s_{\mathrm{xyi},j}^2}},s_{\mathrm{xyi},j}^{n+1},=s_{\mathrm{xyi},j}^{n+1}\frac{\sqrt{\frac{2}{3}}{Y}_0}{\sqrt{s_{\mathrm{xi},j}^2+s_{\mathrm{yi},j}^2+s_{\mathrm{xyi},j}^2}}$

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

附图说明

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

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

具体实施方式

为了说明本发明的具体实施方式，下面将对一个算例进行演示。考虑在二维平面中的问题，一个半无限长的铝棒冲击一个半无限大的铝块，其中左侧铝棒的无量纲初值为u_L＝20.0，p_L＝0.0,ρ_L＝2.7，s_xL＝0.0,s_yL＝0.0，s_xyL＝0.0。右侧铝块的无量纲初值为u_R＝0.0,p_R＝1.0,ρ_R＝2.7,s_xR＝0.0,s_yR＝0.0,s_xyR＝0.0。该问题的无量纲的求解区域为其中铝棒的初始区域为x×y∈[-0.02,0.00]×[-0.006,0.006]，铝块的初始区域为x×y∈[0.00,0.02]×[-0.02,0.02]。同时，铝的理想弹塑性模型的相关无量纲参数分别为ρ₀＝2.71，c₀＝538.0，γ_s＝2.71，K＝740000.0，μ＝265000.0,Y₀＝3000.0。

该问题将在界面左右两侧同时产生弹性波和塑性波。取时间步长Δt＝0.0000001，x方向和y方向的空间步长均为0.00001。采用二维Lax-Friedrich格式计算，得到在时间分别为0.00005，0.0001，0.00015时铝中的负的x方向总应力如图2至图4所示。

Claims (3)

1.一种在欧拉坐标系下计算二维理想弹塑性固体的技术，其特征在于，该技术是一套完整的在欧拉坐标系下计算二维理想弹塑性固体相关物理量(密度速度压力偏应力 )的技术。

2.如权利要求1所述的二维理想弹塑性固体，其特征在于，当其处于弹性状态，满足Hooke定律

$\stackrel{\&CenterDot;}{\mathrm{\&rho;}}=K\frac{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{\mathrm{\&rho;}}$

${\stackrel{\&CenterDot;}{s}}_x=2\mathrm{\&mu;}({\stackrel{\&CenterDot;}{s}}_x+\frac{1}{3}\frac{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{\mathrm{\&rho;}}),{\stackrel{\&CenterDot;}{s}}_y=\mathrm{2\&mu;}({\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_y+\frac{1}{3}\frac{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{\mathrm{\&rho;}}),{\stackrel{\&CenterDot;}{s}}_{\mathrm{xy}}=\mathrm{\&mu;}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{xy}}$

其中K是体积模量，μ是剪切模量，ε_x、ε_y、ε_xy分别是相对应的应变；当其处于塑性状态，满足关系

$p=c_0^2(\mathrm{\&rho;}-{\mathrm{\&rho;}}_0)+({\mathrm{\&gamma;}}_s-1)\mathrm{\&rho;e}$

$s_x=s_x\frac{\sqrt{\frac{2}{3}}{Y}_0}{\sqrt{{s_x}^2+{s_y}^2+{s_{\mathrm{xy}}}^2}},s_y=s_y\frac{\sqrt{\frac{2}{3}}{Y}_0}{\sqrt{{s_x}^2+{s_y}^2+{s_{\mathrm{xy}}}^2}},s_{\mathrm{xy}}=s_{\mathrm{xy}}\frac{\sqrt{\frac{2}{3}}{Y}_0}{\sqrt{{s_x}^2+{s_y}^2+{s_{\mathrm{xy}}}^2}}$

其中c₀，ρ₀，γ_s均为与具体固体有关的常数，Y₀是屈服强度，并且该固体的屈服条件为von Mises屈服条件，即当其偏应力满足

${s_x}^2+{s_y}^2+2{s_{\mathrm{xy}}}^2\&le;{\left(\frac{2}{3}{Y}_0\right)}^2$

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

3.如权利要求1所述的二维理想弹塑性固体的计算技术，其特征在于，Hooke定律中的导数在欧拉坐标系下需要转换成物质导数计算，具体离散形式表示为

$p_{i,j}^{n+1}=p_{\mathrm{old}}+K\ln \left(\frac{{\mathrm{\&rho;}}_{i,j}^{n+1}}{{\mathrm{\&rho;}}_{\mathrm{old}}}\right)$

及

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

其中，ρ_old、p_old、s_xold分别表示欧拉坐标系下的网格点()在上一时间步位置(x_old,y_old)(这里处的密度值、压力值、偏应力值，它们均通过类似如下形式的迎风线性插值得到

$\left\{\begin{array}{c}{\mathrm{\&rho;}}_{\mathrm{old}}=\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)-\frac{\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)-\mathrm{\&rho;}(x_{\mathrm{old}},y_{j-1}^n)}{\mathrm{\&Delta;y}}(y_j^n-y_{\mathrm{old}}),v_{i,j}^{n+1}\&GreaterEqual;0\\ {\mathrm{\&rho;}}_{\mathrm{old}}=\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)+\frac{\mathrm{\&rho;}(x_{\mathrm{old}},y_{j+1}^n)-\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)}{\mathrm{\&Delta;y}}(y_{\mathrm{old}}-y_j^n),v_{i,j}^{n+1}<0\end{array}\right.$

这里

$\left\{\begin{array}{c}\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)={\mathrm{\&rho;}}_{i,j}^n-\frac{{\mathrm{\&rho;}}_{i,j}^n-{\mathrm{\&rho;}}_{i-1,j}^n}{\mathrm{\&Delta;x}}(x_i-x_{\mathrm{old}})\\ \mathrm{\&rho;}(x_{\mathrm{old}},y_{j-1}^n)={\mathrm{\&rho;}}_{i,j-1}^n-\frac{{\mathrm{\&rho;}}_{i,j-1}^n-{\mathrm{\&rho;}}_{i-1,j-1}^n}{\mathrm{\&Delta;x}}(x_i-x_{\mathrm{old}}),u_{i,j}^{n+1}\&GreaterEqual;0\\ \mathrm{\&rho;}(x_{\mathrm{old}},y_{j+1}^n)={\mathrm{\&rho;}}_{i,j+1}^n-\frac{{\mathrm{\&rho;}}_{i,j+1}^n-{\mathrm{\&rho;}}_{i-1,j+1}^n}{\mathrm{\&Delta;x}}(x_i-x_{\mathrm{old}})\end{array}\right.$

$\left\{\begin{array}{c}\mathrm{\&rho;}(x_{\mathrm{old}},y_j^n)={\mathrm{\&rho;}}_{i,j}^n+\frac{{\mathrm{\&rho;}}_{i+1,j}^n-{\mathrm{\&rho;}}_{i-1,j}^n}{\mathrm{\&Delta;x}}(x_{\mathrm{old}}-x_i)\\ \mathrm{\&rho;}(x_{\mathrm{old}},y_{j-1}^n)={\mathrm{\&rho;}}_{i,j-1}^n+\frac{{\mathrm{\&rho;}}_{i+1,j-1}^n-{\mathrm{\&rho;}}_{i,j-1}^n}{\mathrm{\&Delta;x}}(x_{\mathrm{old}}-x_i),u_{i,j}^{n+1}<0\\ \mathrm{\&rho;}(x_{\mathrm{old}},y_{j+1}^n)={\mathrm{\&rho;}}_{i,j+1}^n+\frac{{\mathrm{\&rho;}}_{i+1,j+1}^n-{\mathrm{\&rho;}}_{i,j+1}^n}{\mathrm{\&Delta;x}}(x_{\mathrm{old}}-x_i)\end{array}\right.$

p_old，s_xold，s_yold，s_xyold也通过类似方式进行计算，同时，可按上述的计算过程计算出 $s_{\mathrm{xyi},j}^{n+1}.$