船舶力学

含初始缺陷冰板破碎的双参数近场动力学模拟

PDF下载

李双 ¹^,²^,⁴ , 吕海宁 ¹^,² , 黄小华 ³ , 杨建民 ¹^,²

船舶力学 | 结构力学 2026,30(1): 78-88

收起

船舶力学 |结构力学 2026 , 30 (1) : 78 -88

含初始缺陷冰板破碎的双参数近场动力学模拟

全屏

李双¹^,²^,⁴, 吕海宁¹^,², 黄小华³, 杨建民¹^,²

作者信息

^1.上海交通大学海南研究院，海南三亚 572000

^2.上海交通大学海洋工程国家重点实验室，上海 200240

^3.广西大学防灾与结构安全教育部重点实验室，南宁 530004

^4.江苏科技大学船舶与海洋工程学院，江苏镇江 212100

通讯作者:

吕海宁（1979–），男，副教授，博士生导师，通讯作者，E-mail：haining@sjtu.edu.cn

作者简介:

李　双（1992–），男，博士研究生

Two-parameter peridynamic simulation of failure in ice plate with initial defects

Shuang LI¹^,²^,⁴, Hai-ning LÜ¹^,², Xiao-hua HUANG³, Jian-min YANG¹^,²

Affiliations

^1.Hainan Research Institute, Shanghai Jiao Tong University, Sanya 572000, China

^2.State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

^3.Key Laboratory of Disaster Prevention and Structural Safety of Ministry of Education, Guangxi University, Nanning 530004, China

^4.School of Naval Architecture and Ocean Engineering, Jiangsu University of Science and Technology, Zhenjiang 212100, China

出版时间: 2026-01-15 doi: 10.3969/j.issn.1007-7294.2026.01.009

文章导航

摘要

收起

海冰在凝结过程中不可避免地会产生裂缝或气泡等初始缺陷，这些缺陷会影响冰的力学性能。为了研究含初始缺陷冰板的裂纹扩展和破碎特性，本文基于近场动力学（Peridynamics, PD）理论，引入键的切向刚度系数，发展了双参数PD方法，建立了冰的双参数PD弹脆性模型。该模型突破了传统键基PD固定泊松比的限制，同时继承了键基 PD 的简单性和稳定性等优点。通过对薄板的连续变形和冰梁的三点弯曲断裂试验进行模拟，验证了模型的准确性。研究了冰的泊松比、初始裂隙缺陷、气泡缺陷、冰晶粒尺寸、初始缺陷形状以及多孔缺陷对冰力学特性的影响。结果表明，采用冰的双参数PD弹脆性模型来预测冰体从连续到非连续的全破坏过程是可行的，且冰体中的初始缺陷极大影响着冰体结构的起裂时刻和破坏模式。

关键词

近场动力学 / 冰 / 初始缺陷 / 裂纹扩展 / 脆性破坏

Abstract

收起

Initial defects such as cracks or bubbles will inevitably occur during sea ice condensation, which affect the mechanical properties of ice. To study the crack propagation and fracture characteristics of ice plate with initial defects, based on the peridynamic (PD) theory, the tangential stiffness coefficient for bond is introduced, a two-parameter PD method is proposed, and the two-parameter PD elastic-brittle model of ice material is established. The new model overcomes the limitation of fixing Poisson's ratio and inherits simplicity and stability of the traditional bond-based PD method. The model is verified by simulating the continuous deformation of thin plate and the three-point bending test of sea ice. Furthermore, effects of Poisson’s ratio, initial crack defects, bubble defects, ice grain size, defect shapes and multi-bubble defects on the mechanical properties of sea ice are investigated. The results indicate that the two-parameter PD elastic-brittle model is feasible to predict the continuous to discontinuous failure of sea ice, and the initial defects greatly affect the initiation time and failure mode of the ice structure.

Key words

peridynamics / ice / initial defect / crack propagation / brittle failure

引用本文

李双, 吕海宁, 黄小华, 杨建民. 含初始缺陷冰板破碎的双参数近场动力学模拟. 船舶力学, 2026 , 30 (1) : 78 -88 . DOI: 10.3969/j.issn.1007-7294.2026.01.009

Shuang LI, Hai-ning LÜ, Xiao-hua HUANG, Jian-min YANG. Two-parameter peridynamic simulation of failure in ice plate with initial defects[J]. Journal of Ship Mechanics, 2026 , 30 (1) : 78 -88 . DOI: 10.3969/j.issn.1007-7294.2026.01.009

正文

收起

0　引言

收起

海洋强国战略的提出，加快了我国对海洋极地资源的探索和开发。在此过程中，无论是极地船舶在航行时还是海工结构在作业时，都不可避免地会与海冰发生相互作用，这其中冰荷载的大小和冰体结构的断裂破碎特性是影响结构安全的关键因素^[1]。因此，研究冰的力学性能和破坏特性是支持极地冰区结构物发展、保障极地科考船和极地海洋平台安全性的重要手段，而且海冰的力学特性研究也是朱英富等^[2]提出的极地船舶核心关键技术问题之一。

针对海冰在荷载作用下从受力变形至破坏的过程，相关学者采用了理论推导和物理试验的方法进行了研究。吴刚等^[3]乘坐“雪龙 2”号破冰船在南极考察站完成了冰弯曲强度试验，获取了中山站附近冰的弯曲强度力学特性。余朝歌等^[4]研究了加载速率以及回温对模型海冰劈裂抗拉强度的影响。刚旭皓等^[5]结合原位悬臂梁方法对盐水柱状模型冰的弯曲强度进行了试验研究，测量了模型冰的密度、温度、盐度等基础特性和不同试验条件下的弯曲强度。

相较于成本高昂的物理试验，数值分析方法具有成本低廉、可重复性好等特点，已成为研究冰体力学行为的有效手段。蔡伟等^[6]选取了三种模型试验工况，开展了不同冰材料模型下有限元数值模拟。李晔^[1]为研究海冰与锥体结构的相互作用过程，建立了适用于模拟海冰破碎特性的离散单元模型。雷建奇等^[7]采用粘结离散单元法，研究了加载方位、冰厚度和加载速率对海冰破坏模式的影响。上述有限元方法在不连续处存在奇异性，计算结果具有网格依赖性；离散元法虽然在模拟破坏时表现出优势，但离散单元间采用梁模型，对连续体阶段的计算往往存在较大误差^[8]。

新兴的近场动力学（Peridynamics, PD）^[9] 理论实现了对传统经典连续介质力学的重构，该方法基于非局部作用思想，采用积分型运动方程代替传统微分方程，避免了在不连续处求解微分方程遇到的奇异性问题。并且PD的断裂准则通过一个记忆函数嵌入控制方程中，使之能够自发地模拟裂纹的萌生、扩展至破坏的全过程。在冰的破坏模拟方面，薛彦卓等^[10]采用键基PD方法，建立了冰材料的键基PD弹脆性本构模型，模拟了冰三点弯曲断裂。Wang等^[11]基于键基PD方法，模拟了冰与加筋板接触变形行为。张媛等^[12]采用态基PD方法，建立了冰弹塑性本构模型，模拟了冰梁的四点弯曲破坏。这些基于键的PD冰模型存在固定泊松比的限制（平面应力时，泊松比固定为1/3，平面应变和三维时，泊松比固定为1/4）。基于态基PD的冰材料模型虽然突破了固定泊松比的限制，但态基PD的多体计算方式仍然具有计算复杂和成本高昂等缺点。Zhu 等^[13]通过考虑键旋转来增加键的转角，突破了固定泊松比的限制。与本文模型相比，两者引入切向刚度的方法和推导过程不相同，且Zhu 等未就改进PD 模型开展破坏分析，本文通过在键内直接增加切向刚度系数，在小变形条件下，以均匀应变工况，得到模型的双参数公式，推导过程更加简便，并基于该模型开展了冰板破坏分析。目前，针对冰的裂纹扩展和破碎特性，尤其是含初始缺陷冰材料（固结过程中以及体积膨胀产生的裂缝气泡等缺陷）破碎特性的PD方法研究仍然较少。

因此，本文基于突破固定泊松比限制的双参数PD方法，建立了冰的双参数PD脆性本构模型，并采用FORTRAN 95计算机语言编写程序，模拟研究了含初始缺陷冰板的裂纹扩展和破坏行为。

1　近场动力学理论

收起

近场动力学^[9]是使用节点力的积分来代替运动方程中的空间导数，物体被离散为一系列带质量的物质点，两物质点之间通过“键”来连接。在初始构型中，物质点 ${\boldsymbol{X}}$

的体积为 ${V_{\boldsymbol{X}}}$

，${\boldsymbol{X'}}$

是物体内部物质点， u 和u′ 分别表示变形之后两物质点对应的位移，则在 $t$

时刻，两物质点相对位置矢量 $ {\boldsymbol{\xi }} $

和相对位移矢量 ${\boldsymbol{\eta }}$

分别表述为 $ {\boldsymbol{\xi }}{\text{ = }}{\boldsymbol{X'}} - {\boldsymbol{X}} $

，${\boldsymbol{\eta }} = {\boldsymbol{u'}} - {\boldsymbol{u}}$

，物质点在变形之后的相对位置用 ${\boldsymbol{y}}' - {\boldsymbol{y}} = {\boldsymbol{\xi }}{\text{ + }}{\boldsymbol{\eta }}$

表示。

在键基PD中，连接键包含了物质点的所有本构信息，则物质点 ${\boldsymbol{X}}$

在 t 时刻的运动方程是

$ \rho \ddot{\boldsymbol{u}}=\int_{H_X} \boldsymbol{f}\left(\boldsymbol{u}^{\prime}-\boldsymbol{u}, \boldsymbol{X}^{\prime}-\boldsymbol{X}, t\right) \mathrm{d} V_{X^{\prime}}+\boldsymbol{b}(\boldsymbol{X}, t) $

(1)

式中：$ \rho $

是物质的质量密度；ü 表示加速度矢量场；f 是连接物质点 ${\boldsymbol{X}}$

和 ${\boldsymbol{X'}}$

的对点力；b 是施加的体力密度；${V_{{\boldsymbol{X'}}}}$

表示物质点 ${\boldsymbol{X'}}$

占有的空间体积；内部子区域 ${H_{\boldsymbol{X}}}$

表示以物质点 ${\boldsymbol{X}}$

为中心，以 $\delta $

为半径的近场范围（各物质点在 $\delta $

内相互作用），又称物质点 ${\boldsymbol{X}}$

的“域”，被定义为 $ {H_{\boldsymbol{X}}} = \left\{ {{\boldsymbol{X'}} \in {{\boldsymbol{R}}_{\boldsymbol{0}}}:\left| {{\boldsymbol{X}} - {\boldsymbol{X'}}} \right| \leqslant \delta } \right\} $

。

对线弹性材料而言，物质点的对点力 f 是键的微观弹性应变能密度 $\omega ({\boldsymbol{\xi }},{\boldsymbol{\eta }})$

对该键相对位移矢量 ${\boldsymbol{\eta }}$

的导数^[9]

$ {\boldsymbol{f}}({\boldsymbol{\eta }},{\boldsymbol{\xi }}) = \frac{{\partial \omega ({\boldsymbol{\eta }},{\boldsymbol{\xi }})}}{{\partial {\boldsymbol{\eta }}}} = \frac{1}{2}\frac{{\partial \left( {c({\boldsymbol{\xi }},\delta ){s^2}\left| {\boldsymbol{\xi }} \right|} \right)}}{{\partial {\boldsymbol{\eta }}}} $

(2)

式中：$c$

为键的微观模量；$\left| {\boldsymbol{\xi }} \right|$

为键的原长；$s$

为 t 时刻键的相对伸长率，$s = \left( {\left| {{\boldsymbol{\xi }} + {\boldsymbol{\eta }}} \right| - \left| {\boldsymbol{\xi }} \right|} \right)/\left| {\boldsymbol{\xi }} \right|$

。

在PD模型中，系统的弹性应变能可以通过对微势函数积分得到，将所有连接物质点 ${\boldsymbol{X}}$

的键的应变能密度等分给键两端的物质点，则物质点 ${\boldsymbol{X}}$

的宏观弹性应变能密度为

$ {W^{{\text{PD}}}}({\boldsymbol{X}}) = \frac{1}{2}\int_{{H_\delta }} {\omega ({\boldsymbol{\xi }},{\boldsymbol{\eta }}){\text{d}}{V_{X'}}} $

(3)

在相同的应力或应变状态下，PD物质点 ${\boldsymbol{X}}$

处的应变能密度 $ {W^{{\text{PD}}}}({\boldsymbol{X}}) $

与传统连续介质力学应变能密度 $ {W^{{\text{CCM}}}}({\boldsymbol{X}}) $

是相等的。由此可推导出材料的宏观参数和PD键的微观模量 $c$

之间的关系

$ \left\{\begin{array}{lll}c=\dfrac{9E}{{\text{π}} {\delta }^{3}}&\nu =\dfrac{1}{3},&平面应力\\ c=\dfrac{48E}{{\text{π}} {\delta }^{3}}&\nu =\dfrac{1}{4}\text{，}&平面应变\end{array} \right.$

(4)

式中：E是材料的弹性模量，$v$

是泊松比。这种传统键基PD具有材料泊松比固定的局限。

2　冰板破碎的双参数近场动力学模型

收起

2.1　冰的脆性特性

海冰在凝结过程中体积膨胀，不可避免地产生初始缺陷如裂缝气泡等，这些初始缺陷都会对冰体的断裂破坏产生影响。冰在受力变形至破坏过程中，受加载速率的影响，会出现韧脆转换现象，当加载速度缓慢时，冰表现为韧性材料；高速加载时（应变率大于 10⁻³ s⁻¹），冰表现为极易碎的弹脆性材料^[14]。在冰工程中，海冰的脆性破坏是最受工程人员关注的问题^[11]，因此，本文以高速率荷载下含初始缺陷冰板为研究对象，建立海冰的线弹性双参数PD模型，模拟海冰的裂纹扩展和脆性破坏力学行为。

2.2　双参数近场动力学

为了突破固定泊松比的限制，基于常规微观弹脆性本构模型，在法向刚度的基础上增加切向刚度，提出双参数PD模型。二维情形下，法向刚度对应平行键的初始方向（${{\boldsymbol{\hat e}}_{\rm{n}}}\parallel {\boldsymbol{\xi }}$

），切向刚度对应垂直于键初始方向（${{\boldsymbol{\hat e}}_{\rm{t}}} \bot {\boldsymbol{\xi }}$

），此时，$ {\boldsymbol{X}} $

和 $ {\boldsymbol{X'}} $

间键的微观弹性应变密度为

$ \omega = \frac{1}{2}{c_{\mathrm{n}}}{\left(\frac{{{{\boldsymbol{\eta }}_{\rm{n}}}}}{{\boldsymbol{\xi }}}\right)^2}\xi + \frac{1}{2}{c_{\mathrm{t}}}{\left(\frac{{{{\boldsymbol{\eta }}_{\rm{t}}}}}{{\boldsymbol{\xi }}}\right)^2}{\boldsymbol{\xi }} $

(5)

式中：${c_{\mathrm{n}}}$

和 ${c_{\mathrm{t}}}$

分别为法向和切向刚度系数；${{\boldsymbol{\eta }}_{\rm{n}}}$

和 ${{\boldsymbol{\eta }}_{\rm{t}}}$

分别是键的相对位移矢量，是 ${\boldsymbol{\eta }}$

沿着键初始方向 ${{\boldsymbol{\hat e}}_{\rm{n}}}$

和垂直于该方向 ${{\boldsymbol{\hat e}}_{\rm{t}}}$

的分量；${{\boldsymbol{\eta }}_{\rm{n}}}/{\boldsymbol{\xi }}$

和 ${{\boldsymbol{\eta }}_{\rm{t}}}/{\boldsymbol{\xi }}$

分别表示键在法向和切向方向的伸长率；${{\boldsymbol{\hat e}}_{\rm{n}}}$

和 ${{\boldsymbol{\hat e}}_{\rm{t}}}$

分别表示沿键初始方向和垂直于该方向的单位矢量。

对点力 ${\boldsymbol{f}}$

是键的微观弹性应变能密度对该键相对位移矢量的导数

$ {\boldsymbol{f}} = \frac{{\partial \omega ({\boldsymbol{\eta }},{\boldsymbol{\xi }})}}{{\partial {\boldsymbol{\eta }}}} = \frac{{{c_{\mathrm{n}}}{{\boldsymbol{\eta }}_{\rm{n}}}{{{\boldsymbol{\hat e}}}_{\rm{n}}} + {c_{\mathrm{t}}}{{\boldsymbol{\eta }}_{\rm{t}}}{{{\boldsymbol{\hat e}}}_{\rm{t}}}}}{{\boldsymbol{\xi }}} $

(6)

假定弹性体在小变形情况下发生均匀应变 ${\varepsilon _{ij}}$

，根据Cauchy-Born准则，位移矢量 $ {\boldsymbol{\eta }} $

按法向 $ {{\boldsymbol{\hat e}}_{\rm{n}}} $

和切向 $ {{\boldsymbol{\hat e}}_{\rm{t}}} $

可分解为 $ {{\boldsymbol{\eta }}_{\rm{n}}} $

和 $ {{\boldsymbol{\eta }}_{\mathrm{t}}} $

，因此物质点 ${\boldsymbol{X}}$

的宏观弹性应变能密度为

$ {{\text{W}}^{{\text{PD}}}}({\boldsymbol{X}}) = \frac{1}{2}\int_{{{H}_{\boldsymbol{X}}}} {\frac{1}{{2\xi }}({c_{\mathrm{n}}}\eta _n^2 + {c_{\mathrm{t}}}\eta _t^2)\,{\text{d}}{V_{{\boldsymbol{X'}}}}} $

(7)

化简后可得

$ {W^{{\text{PD}}}}({\boldsymbol{X}}) = \left( {\frac{{{\text{π}} h{\delta ^3}{c_{\mathrm{n}}}}}{{16}} + \frac{{{\text{π}} h{\delta ^3}{c_{\mathrm{t}}}}}{{48}}} \right)\left( {\varepsilon _{11}^2 + \varepsilon _{22}^2} \right) + \left( {\frac{{{\text{π}} h{\delta ^3}{c_{\mathrm{n}}}}}{{24}} - \frac{{{\text{π}} h{\delta ^3}{c_{\mathrm{t}}}}}{{24}}} \right){\varepsilon _{11}}{\varepsilon _{22}} + \left( {\frac{{{\text{π}} h{\delta ^3}{c_{\mathrm{n}}}}}{{12}} + \frac{{{\text{π}} h{\delta ^3}{c_{\mathrm{t}}}}}{{12}}} \right)\varepsilon _{12}^2 $

(8)

在经典连续介质力学中，平面应力状态和平面应变状态的应变能密度分别为

$ {W^{{\text{CCM}}}}({\boldsymbol{X}}) = \frac{E}{{{\text{2(}}1 - {\nu ^2}{\text{)}}}}\left( {\varepsilon _{11}^2 + \varepsilon _{22}^2} \right) + \frac{{E \nu }}{{1 - {\nu ^2}}}{\varepsilon _{11}}{\varepsilon _{22}} + \frac{E}{{1 + \nu }}\varepsilon _{12}^2 $

(9)

$ {W^{{\text{CCM}}}}(X){\text{ = }}\frac{{E(1 - \nu )}}{{2(1{\text{ + }}\nu )(1 - 2\nu )}}\left( {\varepsilon _{11}^2 + \varepsilon _{22}^2} \right){\text{ + }}\frac{{E\nu }}{{(1{\text{ + }}\nu )(1 - 2\nu )}}{\varepsilon _{11}}{\varepsilon _{22}}{\text{ + }}\frac{E}{{(1{\text{ + }}\nu )}}\varepsilon _{12}^2 $

(10)

PD是基于连续介质力学理论框架建立的物质模型，因此对于相同的应力状态，它和连续介质力学的应变能密度必须相等。因此可得二维情形下双参数PD法向和切向刚度系数

$ \left\{\begin{array}{lll}{c}_{{\mathrm{n}}}=\dfrac{6E}{{\text{π}} h{\delta }^{3}(1-\nu )}& {c}_{{\mathrm{t}}}=\dfrac{6E(1-3\nu )}{{\text{π}} h{\delta }^{3}(1-{\nu }^{2})} \text{，}& 平面应力\\ {c}_{{\mathrm{n}}}=\dfrac{6E}{{\text{π}} h{\delta }^{3}(1+\nu )(1-2\nu )}& {c}_{{\mathrm{t}}}=\dfrac{6E(1-4\nu )}{{\text{π}} h{\delta }^{3}(1+\nu )(1-2\nu )}\text{，}& 平面应变\end{array}\right. $

(11)

当切向刚度系数 ${c_{\mathrm{t}}} = 0$

时，双参数PD的微模量参数退化为传统键基PD参数。微模量参数 ${c_{\mathrm{n}}}$

和 ${c_{\mathrm{t}}}$

作为PD的刚度系数，必须满足 ${c_{\mathrm{n}}} \geqslant 0$

且 ${c_{\mathrm{t}}} \geqslant 0$

。因此双参数PD模型的泊松比适用范围分别是，平面应变时 $0 \leqslant \nu \leqslant {1 \mathord{\left/ {\vphantom {1 4}} \right. } 4}$

，平面应力时 $0 \leqslant \nu \leqslant {1 \mathord{\left/ {\vphantom {1 3}} \right. } 3}$

。

2.3　冰体脆性破坏准则

相关研究者^{[10-11, 15]}已证明采用PD脆性材料模型对冰破坏模拟是合适和有效的。在冰的双参数PD模型中，引入“键的临界伸长率”作为张拉破坏准则。当两个物质点键的拉伸 $s$

超过其临界值 ${s_{\mathrm{t}}}$

时，意味着键发生了无法逆转的断裂，此时物质点的相互作用力消失，且压缩极限 ${s_{\mathrm{c}}}$

可认为是拉伸极限 ${s_{\mathrm{t}}}$

的四倍^[10]。为了能够在本构力函数中体现出材料的破坏情况，引入一个与时间 t 有关的标量函数$ \mu \left( {{\boldsymbol{\xi }},{\text{ }}t} \right) $

，当键断裂时$ \mu \left( {{\boldsymbol{\xi }},{\text{ }}t} \right) $

的值从 1 瞬间转为 0，即该标量函数可表示为

$ \mu \left(\xi ,\text{ }t\right)=\left\{\begin{array}{ll}1&\left|s(\xi ) \lt {s}_{0}\right|\text{, 0}\le t\text{'}\le t\\ 0&其他\text{ }\end{array}\right. $

(12)

式中：键的临界拉伸值 ${s_0}$

可以由临界能量释放率 ${G_0}$

求得，在平面应力状态下，$ {s_0} = \sqrt {4{G_0}/h{\delta ^4}{c_{\mathrm{n}}}} $

，因而本构力函数式（6）可改写为

$ {\boldsymbol{f}} = \frac{{{c_{\mathrm{n}}}{{\boldsymbol{\eta }}_{\rm{n}}}{{{\boldsymbol{\hat e}}}_{\rm{n}}} + {c_{\mathrm{t}}}{{\boldsymbol{\eta }}_{\rm{t}}}{{{\boldsymbol{\hat e}}}_{\rm{t}}}}}{{\boldsymbol{\xi }}}\mu \left( {{\boldsymbol{\xi }},{\text{ }}t} \right) $

(13)

为了表达键断裂的积累，PD理论引入损伤模型 $ \varphi ({\boldsymbol{X}},t) $

表达物质点的断裂程度和宏观裂纹，即某一物质点和其邻域范围内的其他物质点之间，力的消失数占总的相互作用力的比例。损伤模型函数是与时间 t 有关的标量函数，其表达式为

$ \varphi ({\boldsymbol{X}},t) = 1 - \frac{{\mathop \smallint \nolimits_{{H_{\boldsymbol{X}}}} \mu \left( {{{\boldsymbol{X}}^{\boldsymbol{'}}} - {\boldsymbol{X}},t} \right){\text{d}}{V_{{{\boldsymbol{X}}^{\boldsymbol{'}}}}}}}{{\mathop \smallint \nolimits_{{H_{\boldsymbol{X}}}} {\text{d}}{V_{{{\boldsymbol{X}}^{\boldsymbol{'}}}}}}} $

(14)

函数 $ \varphi ({\boldsymbol{X}},t) $

表达了在物质点邻域范围内的损伤程度，其取值范围为 [0,1]。若 $ \varphi ({\boldsymbol{X}},t) = 0 $

则表示键完好，未发生断裂；$ \varphi ({\boldsymbol{X}},t) = 1 $

表示物质点与其邻域范围内的所有物质点连接键均发生断键，中心物质点与其他物质点间的作用力完全失效；当 $ \varphi ({\boldsymbol{X}},t) = 0.5 $

时，裂纹萌发。

为了求解上述模型的PD方程，首先将其离散化，其次采用显式中心差分格式求解动态破坏问题，在该差分格式中，一个数值稳定的时间步长应满足

$ \Delta t \lt \sqrt {2\rho /\left( {\displaystyle\sum_{i = 1}^{{N_{\boldsymbol{X}}}} {c\left( {{\upsilon _{{{\boldsymbol{X}}^{\boldsymbol{'}}}}}{V_{{{\boldsymbol{X}}^{\boldsymbol{'}}}}}} \right)/\left| {{\boldsymbol{X}} - {{\boldsymbol{X}}^{\boldsymbol{'}}}} \right|} } \right)} $

(15)

3　算例验证

收起

3.1　数值解与解析解误差验证

为了验证双参数PD模型计算结果的可靠性和鲁棒性，对该模型进行收敛性分析。采用边长0.5 m的正方形薄板，弹性模量 $ E = 200 $

GPa，密度矢量 $ \rho = 7850 $

kg/m³，泊松比 1/3，将薄板离散为一系列带质量的物质点，邻域半径 $\delta {\text{ = }}3.015\Delta $

，沿矩形短边施加的单轴应力荷载为 $ p = 200 $

MPa。设置物质点间距 $\Delta $

分别等于25 mm、20 mm、10 mm、5 mm、2.5 mm 和2 mm，物质点间距越小代表网格越密，如图1所示，然后计算在不同物质点间距情形下数值计算结果与解析解的最大相对误差，最大相对误差分别为9.99%，8.63%，4.48%，3.35%，1.99%和1.66%。显然，随着网格分隔密度的增加，相对误差逐渐减小，数值解趋近于解析解，因此双参数PD模型具备较好的可靠性和鲁棒性。

为了验证双参数PD模型的计算精度以及突破固定泊松比限制的能力，采用和上例相同的几何和材料模型，分别取泊松比为 1/3 和0.15 进行计算，物质点间距 $\Delta {\text{ = }}10{\text{ mm}}$

，采用键基PD模型和双参数PD模型分别对平板在平面应力条件下的变形进行模拟。图2给出了双参数PD模型和传统键基PD模型计算的泊松比分别为1/3和0.15的薄板沿 x = 0方向的位移分布，以及与解析解结果的对比。由图可知，当泊松比为1/3时，两种模型计算的位移只有微小差异，且与解析解非常吻合；当泊松比等于0.15时，传统PD模型的计算结果开始失真，产生极大误差，而双参数PD计算结果依然与解析解吻合较好。因此，双参数PD可用于计算具有不同泊松比的脆性材料的力学行为。

3.2　冰梁三点弯曲破坏验证

为了验证本文模型模拟海冰脆性破坏的能力，建立冰梁的三点弯曲模型，采用双参数PD进行计算，并将数值结果与Liu等^[16]冰梁三点弯曲模型实验结果进行对比。冰梁的试验和数值模型长度均为700 mm，宽度和高度为75 mm。两个固定支撑点之间的距离为600 mm。在冰模型中间，向下施加0.2 mm/s速度边界条件。采用与Liu等^[16]相同的冰材料参数，冰的密度$\rho = 896{\text{ kg/}}{{\text{m}}^3}$

，弹性模量 $E = 1.8{\text{ GPa}}$

，泊松比为0.25，键的临界拉伸值为 $ {s_{\mathrm{t}}}{\text{ = 6}}{\text{.25}} \times {\text{1}}{{\text{0}}^{ - 4}} $

，模型均匀离散化，物质点间距5 mm。

图3展示了冰梁的y向位移动态分布情况。当t = 200时，此时处于加载初期，位移的最大值出现在梁跨中的加载点处；当 t = 2300 时，跨中的最大位移达到0.22 mm，裂纹开始萌发；当 t = 2800 时，裂纹扩展至梁的中心位置，跨中的最大位移达到0.24 mm；当荷载步达到t = 3000时，冰梁发生断裂，此时梁的跨中位移为 0.32 mm。梁破坏时的跨中挠度值与试验结果吻合较好。图4（a）展示了不同荷载步下冰梁动态破坏的双参数PD模拟结果。与图4（b）所示的试验结果相比，可以发现冰的裂纹萌发位置、扩展路径和断裂模式与试验结果吻合较好，说明本文模型模拟冰破坏问题是可行的。

4　含初始缺陷冰板的裂纹扩展和破碎模拟

收起

冰体在形成或挤压碰撞过程中，其内部会形成各种纹理、杂质和大量的初始缺陷，这些缺陷可影响冰的力学性能。本节研究了冰的初始裂隙缺陷、初始气泡缺陷、冰的结晶大小、冰的初始缺陷形状以及多孔缺陷对冰力学特性的影响。建立含初始缺陷冰板模型时，首先定义冰板的几何模型尺寸为 $0.5{\text{ m}} \times 0.5{\text{ m}} \times {\text{0}}{\text{.005 m}}$

，材料参数和3.2节相同，并将冰板进行均匀离散（物质点间距$\Delta $

），各离散物质点通过连接键相互作用。冰板模型中的初始裂缝缺陷是通过打断一系列物质点之间的连接键形成的，孔洞和气泡初始缺陷是通过删除缺陷位置存在的物质点而形成。在外载作用下，本构方程中的键断裂准则一旦达到破坏条件，损伤开始演化，形成宏观裂缝促使试件断裂破坏。模型邻域半径为 $\delta {\text{ = }}3.015\Delta $

，并在模型的上下两端均匀地施加0.02 m/s 的拉伸速度荷载，采用显式中心差分方案进行求解。

4.1　裂缝缺陷角度和材料泊松比对冰板裂纹扩展的影响

建立含有预制裂隙的冰板模型，裂隙位于冰板中心位置，长度为 $0.1{\text{ m}}$

。裂隙与 $x$

轴水平夹角分别为 0°、30°、45°、60°、和90°。在PD数值模拟中，冰板模型的物质点离散网格尺寸为 $\Delta = 0.005{\text{ m}}$

，总的离散物质点总数为 10000个。图5（a）给出了泊松比为0.25的含裂缝缺陷冰板起裂时刻随预制裂缝角度的变化关系。当预制裂缝角度分别为0°、30°、45°、60°和90°时，对应的起裂时刻分别为 $ 192\text{ μs}、 206\text{ μs}、 222\text{ μs}、 265\text{ μs} $

和 $ 474{\text{ μs}} $

，说明在单轴速度荷载下，冰板预制裂隙与垂直荷载方向的夹角越小，裂纹萌生得越早，冰板越容易破坏。图5（b）给出了含0°裂缝缺陷冰板的起裂时刻随材料泊松比变化的关系。分别取泊松比为0.1、0.2、0.25、0.3 和1/3，计算的冰板起裂时刻分别为 $ 197\text{ μs}、 194\text{ μs}、 192\text{ μs}、 190\text{ μs} $

和 $ 189{\text{ μs}} $

，裂纹扩展路径相似，说明在单轴速度荷载下，含初始缺陷冰板的裂纹起裂时刻随泊松比的增大而提前。

不同初始裂隙角度冰板在最终时刻的断裂破坏情形如图6所示。由图6（a）可知，当初始裂隙角度为0°时，裂纹萌生后沿着初始裂隙进行水平扩展，最终破坏时形成了“一”型贯通裂纹；由图6（b）~（d）可知，当初始裂隙角度为30°至60°时，每条预制裂隙的两个尖端各萌生一条扩展裂纹，最终形成类似于“Z”型贯通裂纹；由图6（e）可知，当初始裂隙角度为90°时，初始缺陷裂隙的两个尖端各产生了两条扩展裂纹，分别向左右边界扩展，形成两条类似于“H”型贯通破坏裂缝。结果表明，不同预制裂隙角度的冰板，断裂模式也不相同，但裂纹的萌生全部始于预制裂隙尖端，并沿着垂直于加载方向进行扩展，直至裂纹到达试件边界，造成冰板断裂破碎。

4.2　气泡缺陷孔径对冰板断裂的影响

建立含有初始气泡缺陷的冰板模型，预制单个圆形气泡缺陷位于冰板中心位置，气泡的孔洞直径分别为 20 mm、30 mm、40 mm、60 mm、70 mm、80 mm、100 mm、140 mm和180 mm。在PD数值模拟中，冰板模型的物质点离散网格尺寸为 $\Delta = 0.005{\text{ m}}$

，离散物质点总数随着孔径的变化而变化。

图7给出了含气泡缺陷冰板的起裂时刻随气泡缺陷孔径变化的关系。当缺陷孔径从 20 mm向70 mm增加时，冰板的裂纹起裂时刻从$ 343{\text{ μs}} $

减少为$ 244{\text{ μs}} $

，冰板的裂纹起裂时刻随着气泡缺陷直径的增加而提前，原因是在孔径小于70 mm时，缺陷孔径的增大，使冰板缺陷面积增大促使承载能力降低，从而导致起裂时间提前；当孔径从70 mm向180 mm递增时，冰板的裂纹起裂时刻从$ 244{\text{ μs}} $

延迟到$ 276{\text{ μs}} $

，裂纹起裂随着缺陷孔径的增加出现延迟的现象，原因是在孔径大于70 mm时，孔洞周围的应力集中度减弱，使得孔径大的薄板起裂略有延迟。经过计算，对于有限边界的薄板，当孔径继续增大时，薄板会迅速损伤并断裂，此时孔洞缺陷导致冰晶粒的大量缺失是削弱冰板承载力的主要原因。因此孔洞缺陷并不总是负面影响，合理的布置和利用孔洞缺陷，可以延迟冰板的裂纹起裂。

部分气泡缺陷孔径冰板的最终破坏情形如图8所示。裂纹的萌生全部起源于孔洞的左右两端，都沿着垂直于加载方向进行扩展，萌生两条水平裂纹，生长至试件边界造成冰板断裂破碎。

4.3　冰晶粒尺度对气泡缺陷冰板破碎的影响

PD方法在求解物质点运动时，是将物体离散为一系列带质量的物质点，因此在考虑不同冰晶粒尺度时，可通过在数值模型中设置不同的物质点大小和间距来实现不同冰晶粒尺度。为了研究小气泡缺陷冰板裂纹扩展的空间尺度影响，在冰板中心预制单个圆形气泡缺陷，半径为 60 mm。选取了五种冰晶粒尺度：$\Delta {\text{ = }}2{\text{ mm}}$

、$\Delta {\text{ = }}2.5{\text{ mm}}$

、$\Delta {\text{ = }}3.4{\text{ mm}}$

、$\Delta {\text{ = }}5{\text{ mm}}$

和 $\Delta {\text{ = }}10{\text{ mm}}$

，对应的晶粒尺度冰板物质点总数分别为2469个、9888个、22244个、39552个和61832个，图9给出了缺陷冰板的起裂时刻随冰晶粒尺度变化的关系。上述五种冰晶粒尺度对应的起裂时刻分别为 $ 608{\text{ μs}} $

、$ 541{\text{ μs}} $

、$ 408{\text{ μs}} $

、$ 247{\text{ μs}} $

和 $ 152{\text{ μs}} $

，冰板的裂纹起裂时刻随着冰晶粒尺寸的增加而显著提前。尤其是晶粒尺寸从$\Delta {{ = }}2{\text{ mm}}$

增大到 $\Delta {\text{ = }}5{\text{ mm}}$

时，冰板的承载能力陡然下降，说明冰晶粒越小的冰板抵抗外载的能力越强。

不同冰晶粒尺寸的冰板在最终时刻的断裂破坏情形如图10所示。可以看出，虽然冰晶粒尺度不同，但是最终的冰板破坏模式几乎相同。冰晶粒越大，冰板整体的晶粒总数越少，裂纹的萌生均起源于气泡缺陷的左右两端，沿着垂直于加载方向扩展，直至试件边界，造成冰板断裂破碎。

4.4　不同初始缺陷形状的冰板裂纹扩展模拟

建立含有不同初始缺陷形状的冰板模型，缺陷位于冰板中心位置，分别为纵向矩形、正方形、圆形、椭圆形和菱形。其中矩形的长边方向与椭圆形的短轴方向都与荷载方向一致。矩形缺陷的宽度、正方形缺陷的边长、圆形缺陷的半径、椭圆形缺陷的短轴半径和菱形缺陷的对角线长度均保持一致，为60 mm。在PD数值模拟中，冰板模型的物质点离散网格尺寸为 $\Delta = 0.005{\text{ m}}$

，五种不同形状缺陷的冰板物质点总数分别为9760个、9856个、9888个、9768个和9872个。

图11展示了不同初始缺陷形状冰板起裂时刻的柱状图。由图可知，内含纵向矩形、正方形、圆形、椭圆形和菱形缺陷的冰板的起裂时刻分别为$ 324\text{ μs}、\text{ }306\text{ μs} $

、$ 247\text{ μs}、\text{ }220\text{ μs} $

和$ 212{\text{ μs}} $

。因此，在单轴冲击荷载下，上述几种初始缺陷冰板产生裂纹的先后顺序为菱形 > 椭圆形 > 圆形 > 正方形 > 纵向矩形。

图12给出了不同初始缺陷形状的冰板最终时刻的破坏模式。含纵向矩形和正方形缺陷的冰板，应力主要集中在其所含缺陷的四个角点处，因此裂纹首先从四个角点开始萌发，形成两条贯通裂缝，矩形和方形缺陷冰板两者的冰裂破坏模式类似；含圆形、椭圆形和菱形缺陷的冰板，其裂纹在圆孔左右端或菱形的左右尖端萌发，再沿水平方向扩展直至破坏，最终形成一条贯通裂缝，这三者的破坏模式也基本一致。

4.5　多气泡缺陷的冰板裂纹扩展模拟

多气泡缺陷是冰体中缺陷存在的常态，在冰板中随机布置6个不同尺寸的气泡缺陷，半径分别为25 mm、30 mm、40 mm、40 mm、60 mm 和80 mm。

图13给出了多个时步下（t = 200、2460、2710和2780）多气泡缺陷冰板的x和y位移随裂纹扩展的动态分布。当t = 200 时，x的最大位移分布在冰板的四个角点，y的最大位移分布在冰板上下加载边界；当t = 2460 时，由于气泡缺陷的存在，x和y的最大位移呈现不均匀分布，x的最大位移主要分布在左右边界，y的最大位移主要分布在上下边界；当t = 2710 时，位移在裂纹处呈现明显的不连续；当t = 2780 时，冰板产生贯通裂缝发生破坏，此时x的最大位移为0.02 mm，y的最大位移为0.11 mm。含有多个气泡缺陷的冰板裂纹动态扩展过程如图14所示。由于气泡缺陷之间的相互影响，裂纹呈现不规则的扩展，相近气泡孔之间的裂纹出现交叉融合，气泡缺陷大小和位置的改变都会使裂纹扩展路径产生差异。

实例分析结果表明了双参数PD模型能够很好地模拟含初始缺陷冰板从整体到破坏的全过程，验证了该方法用于海冰的脆性裂纹扩展和破坏分析是可行的。此外，由于采用空间积分方程求解来模拟裂纹扩展，在双参数PD模拟过程中不存在常规连续介质方法模拟时的网格二次划分或裂尖应力奇异性等问题，随着荷载的施加，裂纹自然起裂、扩展和贯通。

5　结论

收起

本文基于近场动力学理论，采用双参数PD模型模拟了含初始缺陷的冰板裂纹扩展和破坏模式。在传统键基PD 的基础上，引入切向刚度系数，推导了键的法向和切向微模量表达式，提出双参数PD方法；结合海冰的脆性破坏特性，给出了冰的脆性破坏准则和数值计算方法；通过两个数值算例验证了模型的准确性；进一步研究了单轴冲击荷载下冰的泊松比、初始裂隙缺陷、初始气泡缺陷、冰的结晶大小、冰的初始缺陷形状以及冰体中多气泡缺陷对冰力学特性的影响。得出以下结论：

（1）双参数PD键刚度系数在切向刚度为零时，能完全退化为传统键基PD模型参数。双参数PD数值模拟的非破坏变形算例，与解析解对比误差较小；数值模拟的冰梁三点弯曲开裂位置、扩展方向和最终破坏路径与试验结果吻合较好，证明了本文模型预测海冰脆性破坏的适用性和准确性。

（2）初始缺陷对冰板的裂纹萌发时刻和破坏模式产生较大影响，且泊松比的增大会使冰板起裂时间提前。在一定条件下，裂缝缺陷与垂直荷载方向夹角越小，裂纹萌生得越早；冰晶粒尺寸越大，裂纹萌生得越早；冰板裂纹起裂时刻随着气泡缺陷孔径的增大呈现先提前后延迟的现象；几种初始缺陷形状冰板产生裂纹的先后顺序为菱形−椭圆形−圆形−正方形−纵向矩形；对于多孔洞缺陷，裂纹呈现不规则扩展，相近气泡孔之间相互影响，裂纹出现交叉融合的现象，多气泡缺陷大小和分布的稍有改变，都会引起裂纹扩展路径产生差异。

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

李晔. 锥体海洋平台结构冰荷载的离散单元分析[J]. 船舶力学, 2017, 21(10): 1254‒1262.

Li Y. Analysis of ice load on conical structure with discrete element model[J]. Journal of Ship Mechanics, 2017, 21(10): 1254‒1262. (in Chinese)

朱英富, 刘祖源, 解德, 等. 极地船舶核心关键基础技术现状及我国发展对策[J]. 中国科学基金, 2015, 29(3): 178‒186.

Zhu Y F, Liu Z Y, Xie D, et al. Advancements of the core fundamental technologies and strategies of China regarding the research and development on polar ships[J]. Bulletin of National Natural Science Foundation of China, 2015, 29(3): 178‒186. (in Chinese)

吴刚, 唐文勇, 王庆凯, 等. “雪龙2”号南极首航破冰试验[J]. 船舶力学, 2021, 25(8): 981‒990.

Wu G, Tang W Y, Wang Q K, et al. Full-scale ice trials of R/V Xuelong 2 during hermaiden Antarctic voyage[J]. Journal of Ship Mechanics, 2021, 25(8): 981‒990. (in Chinese)

余朝歌, 田于逵, 刚旭皓, 等. 基于巴西圆盘劈裂试验的模型海冰抗拉强度研究[J]. 船舶力学, 2020, 24(9): 1158‒1166.

Yu C G, Tian Y K, Gang X H, et al. Investigation on the tensile strength of model ice based on Brazilian disc splitting test[J]. Journal of Ship Mechanics, 2020, 24(9): 1158‒1166. (in Chinese)

刚旭皓, 田于逵, 季少鹏, 等. 盐水柱状模型冰弯曲强度试验分析[J]. 船舶力学, 2021, 25(3): 336‒341.

Gang X H, Tian Y K, Ji S P, et al. Experimental analysis on flexural strength of columnar saline model ice[J]. Journal of Ship Mechanics, 2021, 25(3): 336‒341. (in Chinese)

蔡伟, 朱凌, 毕璐泽. 船-冰碰撞下冰材料有限元数值方法研究进展[J]. 船舶力学, 2021, 25(8): 1116‒1126.

Cai W, Zhu L, Bi L Z. Review on development of finite element numerical method of ice material model under ship-ice collision[J]. Journal of Ship Mechanics, 2021, 25(8): 1116‒1126. (in Chinese)

雷建奇, 狄少丞, 于加一, 等. 船-冰作用中海冰典型破坏模式的离散元数值模拟[J]. 计算力学学报, 2021, 38(3): 346‒354.

Lei J Q, Di S C, Yu J Y, et al. Discrete element numerical simulation of typical failure modes of sea ice under ship - ice interaction[J]. Chinese Journal of Computational Mechanics, 2021, 38(3): 346‒354. (in Chinese)

Selvadurai A, Sepehr K. Two-dimensional discrete element simulations ofice–structure interaction[J]. International Journal of Solids and Structures, 1999, 36(31-32): 4919‒4940.

Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces[J]. Journal of the Mechanics & Physics of Solids, 2000, 48(1): 175‒209.

薛彦卓, 陆锡奎, 王庆, 等. 冰三点弯曲试验的近场动力学数值模拟[J]. 哈尔滨工程大学学报, 2018, 39(4): 607‒613.

Xu Y Z, Lu X K, Wang Q, et al. Simulation of three-point bending test of ice based on peridynamic[J]. Journal of Harbin Engineering University, 2018, 39(4): 607‒613. (in Chinese)

Wang Q, Wang Y, Yuan L, et al. Simulation of brittle-ice contacting with stiffened plate with peridynamics[J]. Journal of Ship Mechanics, 2018, 22(3): 339‒352.

张媛, 王超, 郭春雨, 等. 冰弯曲破坏的弹塑性近场动力学模型[J]. 哈尔滨工业大学学报, 2022, 54(6): 87‒94.

Zhang Y, Wang C, Guo C Y, et al. An elastoplastic model of ice bending failure in peridynamics[J]. Journal of Harbin Institute of Technology, 2022, 54(6): 87‒94. (in Chinese)

Zhu Q, Ni T. Peridynamic formulations enriched with bond rotation effects[J]. International Journal of Engineering Science, 2017, 121: 118‒129.

Schulson E M. Brittle failure of ice[J]. Engineering Fracture Mechanics, 2001, 68(17-18): 1839‒1887.

Kan X, Zhang A, Yan J, et al. Numerical investigation of ice breaking by a high-pressure bubble based on a coupled BEM-PD model[J]. Journal of Fluids and Structures, 2020, 96: 103016.

Liu R W, Xue Y Z, Lu X K, et al. Simulation of ship navigation in ice rubble based on peridynamics[J]. Ocean Engineering, 2018, 148: 286‒298.

2026年第30卷第1期

PDF下载

112

引用本文

BibTeX

文章信息

doi: 10.3969/j.issn.1007-7294.2026.01.009

接收时间：2025-02-27
首发时间：2026-07-07
出版时间：2026-01-15

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-02-27

基金

作者信息

^1.上海交通大学海南研究院，海南三亚 572000

^2.上海交通大学海洋工程国家重点实验室，上海 200240

^3.广西大学防灾与结构安全教育部重点实验室，南宁 530004

^4.江苏科技大学船舶与海洋工程学院，江苏镇江 212100

通讯作者:

吕海宁（1979–），男，副教授，博士生导师，通讯作者，E-mail：haining@sjtu.edu.cn

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/cblx/CN/10.3969/j.issn.1007-7294.2026.01.009

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT