数值仿真中单元密度对材料失效应变的影响

,

(上海交通大学 船舶海洋与建筑工程学院，上海 200030)

近几年来数值方法在船舶碰撞与搁浅、舰船水下爆炸等方面得到了广泛应用[1-3]。在这些非线性瞬态动力学问题的数值仿真计算中不得不考虑结构塑性失效问题，数值计算中一般通过设定单元最大塑性失效应变εmax来模拟材料失效[4-5]，因此材料最大塑性失效应变εmax的正确选择直接关系到计算结果的准确性。在以往的研究中，材料失效应变如何设定并没有统一规定，有些直接将材料拉伸实验数据用于有限元计算，有些则依据计算经验选取。本文以普通船用低碳钢为例，对材料拉伸实验进行数值仿真校准分析。

1 理论背景

材料应力-应变曲线通常根据试样拉伸实验得到。应力和应变的计算方法如下：

(1)

式中：F——试样拉伸过程中仪器记录的载荷；

A0——试样原始横截面积；

L0——试样原始长度。

根据式(1)计算得到的应力、应变是在试样特征长度(即试样长度)内的平均应力和平均应变，叫工程应力σE和工程应变εE[6]。

随着试样不断被拉伸，试样发生颈缩，随后的拉伸变形主要集中在颈缩区域，颈缩处横截面变得越来越小，相应的局部应变和局部真实应力不断增加，直到试样被拉断。图1所示为应力应变曲线[7]。

图1 应力-应变曲线

工程应力σE、应变εE与真实应力σT、应变εT的关系如下[8]：

σT=σE(1+εE)

(2)

εT=ln(1+εE)

(3)

上述真实应力-应变关系在发生颈缩后不再适用。发生颈缩后材料的真实应力-应变关系采用Hollomon乘幂关系描述[9]：

(4)

数值仿真计算中单元的受力状态见图2。

图2 单元受力状

单元应力、应变为：

(5)

即数值计算中应力、应变的真实含义是：单元所在位置单元内的平均应力和应变。因此当单元长度与拉伸实验试样的特征长度不相等时，不能直接使用拉伸实验得到的断裂应变作为数值计算中材料的最大失效应变εmax。

2 拉伸实验数值仿真校准

在数值仿真中，采用真实应力-应变曲线定义材料的力学性能，将拉伸试样一端固定，一端施加强迫位移，输出材料固定端的固端力，该载荷等效于试样拉伸过程中实验仪器记录的载荷F，根据式(1)即可得到数值仿真计算得到的工程应力、应变，记做σEN和εEN。对比数值仿真计算得到的工程应力-应变曲线和实验得到的应力-应变曲线即可对材料拉伸实验进行数值校准。

数值校准采用的程序为MSC.DYTRAN，令材料拉伸实验得到的工程断裂应变为εEF，数值仿真中材料的最大塑性失效应变为εmax，εmax即为对应于单元特征长度的断裂应变。两者之间的关系用表示如下：

εmax=αεEF

(6)

图3所示为数值仿真校准计算中材料的应力-应变曲线。

图3 试样应力-应变曲

试样几何尺寸如图4所示。数值计算中为了保证材料变形处于准静态范围内，端部强迫位移速度为5 mm/s。共对8种不同网格密度的试样模型进行了数值校准计算。具体网格尺寸参数见表1。有限元模型见图5。

图4 试样几何尺

表1 有限元网格参数

图5 颈缩现象的形成

3 仿真结果分析

3.1 颈缩现象

如前所述，随着试样被不断拉伸，在试样中将出现颈缩想象，颈缩形成前，整个试样内变形是均匀的，因此试样内的应变分布也是均匀的。发生颈缩后变形主要集中在颈缩区域，则颈缩区域的应变比周围区域要大。这即是需要对单元最大失效应变进行校准的原因。

图5给出了Case-1和Case-8两种模型的拉伸变形情况。

Case-1的单元个数为1，其大小与试样几何尺寸相同，在数值仿真拉伸过程中单元应变即为试样的平均应变，拉伸时单元截面因变形而收缩，但不会出现颈缩现象，因此当单元几何尺寸大于或等于拉伸试样特征长度时单元的失效应变即为拉伸实验得到的工程断裂应变，修正系数α≈1。

Case-8的单元大小为0.5 mm，从图5可以看出，当拉伸量ΔL=20 mm时，试样出现了颈缩现象。

图6给出了颈缩后试样内的应变分布，可见在颈缩区域应变分布极不均匀。

图6 应变分布云图

图7为试样中心点处的单元应变随拉伸距离的变化关系。

图7 试样中心点处应变-拉伸量曲

颈缩前不同模型的应变在整个试样范围内均匀分布，约等于ΔL/L。当拉伸量达到ΔL=15 mm后，试样开始出现颈缩想象，由于颈缩的形成，不同模型的曲线出现了急剧发散。由于单元尺寸不同，单元覆盖的区域也不同，每个单元表示的是该区域内的应变大小，即单元越小越接近于该点的真实应变。

3.2 最大塑性失效应变修正

本文试样拉伸实验测得的工程断裂应变为0.368。

表2给出了拉伸数值校准得到的不同网格密度情况下，单元的最大失效应变。

表2 最大失效应变修正

图8为修正曲线。随着单元密度的增加，单元的最大失效应变急剧增加，单元尺寸接近试样特征尺寸时，修正系数趋近于1。

图8 最大失效应变修正系

普通低碳钢的临界破坏应力为1 533 MPa[8-9]，该临界应力即为低碳钢断裂时材料分子间能够承受的最大载荷，由此可见，当单元尺寸减小时，单元的失效应力将逐渐增加；如果单元特征尺寸足够小，单元失效应力将接近材料的临界破坏应力。与此相对应的失效应变将随着单元特征尺寸减小而增加。根据本文计算，当材料单元尺寸达到0.5 mm时，单元失效应变为3.32，对应的材料真实应力为764.5 MPa。

4 结论

1) 由于数值计算中单元应力和应变即为单元所在区域的平均应力和应变，因此，当单元尺寸大于试样特征尺寸时，单元最大失效应变即为拉伸实验测得的工程断裂应变，此时的修正系数α≈1。

2) 当单元尺寸小于试样特征尺寸时，单元应力和应变更接近于单元所在位置的真实应力和真实应变，此时需要根据本文的修正曲线对数值计算中单元最大失效应变进行修正。

[1] 王自力.船舶碰撞损伤机理与结构耐撞性研究[D].上海: 上海交通大学, 2001.

[2] Rajendran R, Narasimhan K.Damage prediction of clamped circular plates subjected to contact underwater explosion [J].International Journal of Impact Engineering, 2001，25:373-386.

[3] 陈永念，尹 群，胡海岩.水中爆炸冲击波载荷作用下舰船结构动态响应的数值模拟 [J].爆炸与冲击,2004，24：201-206.

[4] R Rajendran and K Narasimhan.Deformation and fracture behaviour of plate specimens subjected to underwater explosion-a review [J].International Journal of Impact Engineering， 2006,32:1945-1963.

[5] 梅志远，朱 锡，刘润泉.船用加筋板架爆炸载荷下动态响应数值分析[J].爆炸与冲击, 2004,24:80-84.

[6] Courtney.Mechanical Behavior of Materials [M].New York: McGraw-Hill, 1990.

[7] Boyer.Atlas of Stress-Strain Cures: ASM International, 1987.

[8] Young W C.Roark’s Formulas for Stress and Strain[M].7th ed.McGraw-Hill Companies Inc, 2002.

[9] 盖京波.舰船结构在爆炸冲击载荷作用下的局部破坏研究[D].哈尔滨: 哈尔滨工程大学, 2005.