冰弯曲破坏的弹塑性近场动力学模型

2022-07-02郭春雨叶礼裕

哈尔滨工业大学学报 2022年6期

张媛，王超，郭春雨，叶礼裕，刘正

(哈尔滨工程大学船舶工程学院，哈尔滨 150001)

近些年来，在气候变暖、科学技术发展、经济效益驱动以及一些政治和军事领域的显著价值需求下，极地航线开发、极地资源的勘探以及极地环境利用成为一种必然且迫切的趋势[1]。这些变化极大地促进了更加频繁的极地地区海洋结构物的作业活动，除此之外，多数内陆结冰水域仍旧常年进行各种海上作业活动。然而冰的存在能够对结构物造成严重的损害，甚至造成结构物的破坏而无法正常作业。由此研究冰力学性能和破坏模式及其与结构物相互作用过程是支持极地以及冰区结构物发展的重要手段之一，而建立精确的冰本构模型是这其中的关键技术之一。

基于数值计算的冰力学特性、冰-结构物作用过程等研究近几年快速发展：最初，学者们通常采用基于试验获得的经验数据和数学物理分析的方法来计算冰的力学特性并且用于冰与结构物相互作用过程的冰载荷预报以及船舶运动性能预报[2-4]。然而上述文献中对于海冰性质研究的方法，极大地依赖于经验数据，不足以发展成完整成熟的数值模型[5]。后来，学者们纷纷采用基于传统连续介质力学方法的数值手段来研究海冰的特性，例如有限元方法FEM[6-8]，改进有限元方法[9-10]，以及其他考虑冰细微参数的方法等[11]。然而传统介质力学方法在处理大变形或者断裂问题具有一定的局限性，不能准确地模拟冰破坏过程的裂纹扩展路径。为了克服上述的科学问题，学者开始尝试建立基于无网格方法的冰本构模型，这些方法里发展比较成型的有：离散元方法(Discrete element method,DEM)[12-13]，光滑粒子动力学(Smoothed particle hydrodynamics,SPH)[14-16]以及近场动力学方法(Peridynamics,PD)[17-19]等。其中DEM方法由于其颗粒状的模型形式具有模拟浮碎冰和结构物碰撞、浮碎冰之间碰撞的优势。而SPH和PD方法的粒子离散形式使其在海冰破碎以及裂纹扩展方面更具有更突出的优势，不管是哪种粒子数值模型，虽然能够处理断裂问题，但是并不能完整精确地模拟冰力学性能和破坏特性，由此这些方法仍然处于发展前期，还需要投入大量的研究工作来建立适用于冰的准确数值模型。

近场动力学方法(PD)，是近二十年新兴的一种非局部、无网格粒子、积分形式的计算方法，在断裂问题域和非断裂问题域上都可以连续求解[20]。该方法最早由Silling等[21-22]提出并建立推导了其控制方程，给出了基本的数值实现程序案例。PD方法主要两种表现形式[21-23]，经推导和比较[24]两种表达方式大同小异且具有相类似的计算精度。PD发展至今主要有键型、常规状态型和非常规状态型3种形式。且该方法已经成功应用于材料的弹性、塑性、热力学等多种物理问题计算分析领域[25-27]。然而PD应用于冰计算领域刚刚起步，多数文献采用基于键型的弹性模型来模拟冰的破坏特性，例如，冰与刚体柱的作用过程、冰桨作用过程等，忽略了泊松比的真实状态和海冰具有塑性屈服过程的力学特性。由此本文主要的目的是建立基于状态型PD方法的冰弹塑性本构模型，消除了键型对泊松比的限制，加入了冰破坏过程塑性变形过程，为后续的冰与结构物的作用过程提供支撑。

1 基于常规近场动力学的弹塑性本构模型

1.1 基于常规近场动力学的弹塑性理论

1.1.1 近场动力学理论

近场动力学方法将连续介质离散为均匀的物质点。在参考坐标系下，其运动方程为

t′(u-u′,x′-x,))dH+b(x,t)

(1)

上式的离散形式为

t′(u(k)-u(j),x(k)-x(j),t))dV(j)…+b(k)

(2)

式中，每个粒子的信息用其位置坐标表示，即x(k)，该粒子占据的空间大小为V(k)且密度为ρ(k)。在笛卡尔坐标中，材料点受到外界作用的位移为u(k)，此时的位置向量为y(k)。由此，材料点k在未变形状态下的位置为x(k)，其位移和外力向量分别用u(k)(x(k),t)和b(k)(x(k),t)表示，对k点有作用的域(近场域)表示为Hx(k)。

1.1.2 塑性理论

本文将海冰视作各向同性的均匀介质材料。下述本构模型的建立基于该假设成立。根据经典连续介质力学中的塑性定义，当材料在加载/卸载条件下的应力水平超过屈服应力时，材料在完全卸载时会发生永久(或塑性)变形。与弹性变形不同，塑性变形取决于加载历史。因此，塑性变形更加关注每一时间步力密度增量和拉伸增量的关系。用PD进行塑性分析时，需要对等效应力进行描述，这是因为在塑性分析中，弹性区域受屈服面限制。此外，屈服面在每个加载步骤中的增长和移动都需要确定，因此需要将等效塑性拉伸描述为硬化内变量。

为了用参量的增量形式表达PD中的塑性变形，将外界作用下两个材料点之间的拉伸增量Δs(k)(j)以及力密度增量Δt(k)(j)分解成弹性项和塑性项如下:

(3)

(4)

本文的塑性理论中假设塑性变形与体积变形无关，由此塑性拉伸的过程体积膨胀为0，即Δθp=0。PD方法中静水压力分别在力密度函数和应变能密度的弹性项中，用p=-kθ表述，其中k为冰的体积模量。本文模型中塑性变形与静水压力无关。同样地，每一步的应变能密度增量分解为弹性应变能密度增量和塑性应变能密度增量如下:

(5)

1.1.3 屈服函数

传统的塑性理论中，塑性变形由屈服面来判定。参照传统塑性理论屈服面定义为[28]

(6)

(7)

依据Von-Mises 塑性屈服准则，塑性屈服发生在偏斜应变部分的应变能密度达到极限值的时候，对于单轴压缩来说，这个极限值等于材料达到初始屈服应力σ0时的应变能密度，即

(8)

式中J2为应力偏张量的第二不变量。为了在一般荷载情况下使用该方程，需要定义等效应力。等效应力可定义为[29]

(9)

1.1.4 流动法则和硬化规律

(10)

式中，C(k)为一个正比例常数，上式作为常态条件表述了塑性伸长增量的流动法则。

在式(7)中，理想塑性(无硬化)状态下，所有加载/卸载过程中，σy是恒定的。然而，塑性变形在不同的状态下与等效塑性拉伸具有相关性，这种相关性有多种情况，例如各向同性硬化关系、运动硬化和混合硬化关系等[28]。本文采用线性各向同性硬化模型[29]，即式(6)可改写为

(11)

其中

在简单的单轴拉伸情况下，无膨胀项，可以确定塑性增量拉伸引起的畸变应变能密度，并由此提取参数A0。参数A0已在文献[30]中针对不同的物体维度推导得出，本文未列出。

1.2 边界条件

由于PD运动方程不包含任何空间导数，因此一般不需要约束条件来求解积分微分方程。与局部理论不同，边界条件是通过非零体积的虚拟边界层施加的。基于数值实验，Macek等[30]建议虚拟边界层的范围等于领域尺寸δ，以确保施加的规定约束在真实区域中得到准确反映。位移边界条件可以通过对虚拟层中的材料点进行约束来施加，使边界曲面上的条件得到明确满足。因此，虚拟层中的位移值可以基于实域值和边界条件指定值的线性外推来近似。类似地，也可以通过近似虚拟区域中的牵引力值来施加牵引力边界条件，从而使真实区域和虚拟层中的牵引力变化恢复边界表面上施加的牵引力。图1表示了位移边界条件的施加方式，图中Rf表示虚拟粒子边界，宽度为δ，R为真实粒子区域，位移和速度边界条件为U*和V*。则边界条件可以表达为：

图1 位移和速度约束条件[31]Fig.1 Boundary condition of displacement and velocity[31]

uf(xf,yf,t+Δt)=2U*(x*,y*,t+Δt)-u(x,y,t)

vf(xf,yf,t+Δt)=2V*(x*,y*,t+Δt)-v(x,y,t)

(12)

1.3 失效准则

基于PD方法的弹性变形中，失效的判断标准为极限伸长量，当变形后粒子间的伸长量达到极限伸长量s0时，则判断为失效，此过程不可逆。极限伸长量可以通过材料的极限能量释放率得到，粒子间的伸长量和力密度的关系呈线性。然而在塑性变形中，粒子间力密度和伸长量的关系呈典型的非线性，如图2所示，因此弹塑性本构模型中的失效准则不能再简单地使用最终状态的伸长量来判断，必须对整个变形过程伸长量对应的力密度积分来判断材料的失效。图2曲线下的面积代表弹塑性变形微势w(k)(j)，该函数代表变形的程度，可由下式得到:

图2 塑性理论中的变形本构关系Fig.2 Constitutive relationship in plastic deformation

(13)

如果一个新的裂纹产生，则穿越这个裂纹的所有粒子间的作用都必须消失。同理，消除两个粒子间作用所需的极限能量释放率作为判别粒子间作用是否失效的标准，则可以表达为

(14)

(15)

(16)

式中：Gc为材料的极限能量释放率，Nc为PD的对应参数，在一维时取值为12，二维时取值为36，三维时取值为560。

粒子作用的失效通过上述准则判断之后，失效的粒子作用之间力将被移除。此时，通过监测粒子之间相互作用的失效与否就可以表述材料的破坏过程。因此引入局部破坏系数φ(x(k),t)来表述材料在局部位置的破坏程度，该参数的物理意义是：材料点k在其作用域内失效粒子作用与总粒子作用的比值，由下式计算：

(17)

2 数值执行过程

2.1 数值执行策略

PD方法的控制方程主要采取积分形式，因此在数值实现时，计算域被分解成粒子形式。计算每个粒子在每个时间步长的物理信息，作为下一时间步的基础值，每个时间步下对所有粒子的物理信息求和积分则为该步长的变形状态。

本文参考了文献[24]中给出的键型PD的程序基础框架。本文的数值执行过程在上述程序框架进行补充完善。首先，常规型PD积分过程中增加了体积膨胀项，可按下式数值执行:

(18)

时间积分内的数值执行过程将全部替换为塑性变形的数值过程，可参考文献[28-29]中的数值实现方法。对于失效准则的积分过程，式(13)将采用如下积分过程实现:

t+Δtw(k)(j)=tw(k)(j)+Δw(k)(j)

(19)

(20)

局部破坏系数的数值执行形式为

(21)

数值框架流程如图3所示。

图3 数值实现程序框图Fig.3 Flowchart of numerical implementation

2.2 数值方法验证

为了验证文章模型的准确性和程序的可靠性，首先计算了带孔二维板的变形问题。二维板的尺寸选取长和宽分别为L=H=1.0 m，板的厚度为h=0.01 m，孔的尺寸为D=0.3 m。弹性模量为E=200 MPa，密度为ρ=4 428 kg/m3，泊松比为ν=1/3，刚度模量为K=50 GPa。离散状态的粒子间距为dx=0.002 5 m。同时采用有限元分析软件Abaqus计算了板的有限元变形，有限元分析中，网格大小与粒子间距相同，网格数有8 910个。板的左右两侧采用位移约束，其位移变形约束表示为：ux(x=0,y,t)=-0.001 m和uy(x=L,y,t)=0.001 m。板的工况示意图如图4所示。

图4 二维板拉伸变形示意图Fig.4 Sketch map of 2D plate under tension

正如预期的那样，塑性变形在高应力集中区域开始。图5中给出了带孔板的等效应力和位移的变化，与有限元结果对比完全一致，该案例验证了本文程序的准确性。

图5 有孔板拉伸状态下的等效应力、x方向位移以及y方向位移Fig.5 Variations of effective stress, displacement in x direction, and displacement in y direction of plate with a hole under tension

3 冰四点弯曲数值计算

冰的弯曲强度是计算冰-船作用的关键参数之一，四点弯曲试验可以用来测量冰的弯曲强度。为了验证本文冰模型的准确性，选用Ehlers和Kujala进行的试验[32]作为对比验证，该试验模型及加载工况示意如图6所示。

图6 冰四点弯曲试验示意图Fig.6 Schematic diagram of four-point bending test of ice

该试验中，冰梁的尺寸为L=4.25 m,H=0.4 m和B=0.365 m。上方的固定支撑位于中点向两边2 m处，底部位移压头位于0.5 m处。上方的支撑在空间中固定，下方的支撑以Vsupport=0.003 580 m/s匀速向上移动。海冰参数的选取参考试验和文献中对冰力学性能相关分析设定[32-34]，冰材料的极限能量释放率采用了文献[35]中的参考值，见表1。

表1 试验工况设定和海冰参数Tab.1 Test conditions and parameters of ice

文中采用2D进行计算，因此弹性模量转换成2D参数为E/(1-v2)[34]。结果如图7、8所示。冰的弯曲载荷和时间的关系曲线如图9所示。

从计算结果可以发现，冰梁在下部位移支撑和上层约束支撑的共同作用，发生微小变形，如图7(a)、(b)所示，两个移动支撑之间的冰梁位移最大，且这段部分各个位置的位移量基本相同。冰梁上越远离移动支撑部位的纵向位移变形越小，这样的变形一直持续到0.31 s之前。此时海冰未发生破坏，由于固定支撑和加载支撑的约束原因，冰x方向和y方向的位移符合物理规律。还可以通过等效应力分布图观测到：高应力分布在拉伸和压缩集中的部位，例如冰梁的上侧边缘和下侧边缘，如图5(c)所示，这两侧的位置是冰梁挤压和拉伸最为严重的位置。该应力现象与文献[35]中通过有限元方法计算的结果具有很高的一致性。从图8可以观测到：冰梁弯曲后分裂成3部分，断裂位置对应于底部压头加载位置。该现象和试验现象具有极高的一致性。另外为了进一步验证计算结果的准确性，对比数值计算和试验中的力-时间曲线(图9)可知：数值计算结果和试验结果吻合度极高。

图7 t=0.31 s,Vsupport=0.003 580 m/s(未发生破坏之前)冰四点弯曲试验数值计算结果Fig.7 Numerical results of four-point bending test of ice (before damage) (t=0.31 s,Vsupport=0.003 580 m/s)

图8 t=1.00 s,Vsupport=0.003 580 m/s(破坏之后)冰四点弯曲试验数值计算结果Fig.8 Numerical results of four-point bending test of ice (after damage) (t=1.00 s,Vsupport=0.003 580 m/s)

图9 Vsupport=0.003 580 m/s冰四点弯曲载荷-时间曲线图Fig.9 Force-time curve of four-point bending test of ice (Vsupport=0.003 580 m/s)

为了验证本文的数值模型能够模拟不同时间工况下冰梁的弯曲破坏，分别选取了Vsupport=0.002 750 m/s和Vsupport=0.003 147 m/s加载工况进行对比验证，载荷-时间变化关系计算结果如图10所示，数值结果和试验结果基本保持一致。由试验结果可知：加载速度为Vsupport=0.003 580 m/s时的弯曲最大载荷最大，断裂时间最短。表2列出了3种加载速度下数值计算和试验对应的极值载荷数值极其发生时间，从表中可以看出试验的结果和计算值的误差在10%以内，说明本文建立的无网格法冰本构模型能够模拟冰的弯曲破坏过程。