盐水柱状模型冰弯曲强度数值计算分析

刚旭皓，田于逵，余朝歌，季少鹏，寇莹

中国船舶科学研究中心，江苏 无锡 214082

0 引 言

船舶在冰区中航行时会受到诸多因素的影响，如海冰物理力学特性、海冰密集度、船体形状、航行速度等，其中冰的物理力学特性是通过直接影响冰层的承载能力来对船舶冰阻力产生最本质的影响。在冰区船舶破冰过程中，弯曲破坏为冰层主要的破坏形式之一[1]，因此，研究冰的弯曲强度对于冰区航行船舶的破冰阻力具有重要意义。

自20 世纪70 年代以来，研究人员针对冰的弯曲强度陆续开展了研究，当前针对冰弯曲强度的研究方式主要有现场试验、模型试验和数值计算等方法。Maattanen[2]在波罗的海的现场试验中就冰样尺寸和加载速率对海冰弯曲强度的影响进行了研究，得到冰的弯曲强度会随冰样宽度的增加而增加，但几乎不受冰样长度影响的结论。Frederking和Timco[3-4]分别在不同的时间就模型冰的弯曲强度进行了试验，测得冰的弯曲强度与加载速率无明显相关性。李志军等[5]对细粒酒精冰弯曲强度的影响因素进行研究，确定了加载速率对模型冰弯曲强度的影响与渤海海冰类似，均呈韧脆转变的特性。在数值计算方面，季顺迎[6]利用离散元方法（discrete element method，DEM）对冰的压缩强度和弯曲强度进行了数值计算，并对计算模型的细观参数予以了校准；韩端峰等[7]利用SPH（smoothed particle hydrodynamics）方法对冰的三点弯曲进行模拟，定义了粒子间的作用类型。DEM 法与SPH法均是采用粒子法来对冰层进行模拟，粒子法的关键是确定粒子之间接触的搜索方法。魏龙海等[8]采用链式结构搜索法，利用全局搜索的方法完整记录了链式结构体内的数据。刘璐[9]开发了邻居列表法，在搜索中只需对相邻颗粒进行搜索判断即可，节省了计算时间。

本文将在前期中国船舶科学研究中心小型冰水池研究的基础上，利用原位悬臂梁方法对盐水柱状模型冰的弯曲强度进行一系列试验研究，分析弯曲强度的影响因素，然后在此基础上利用冰力学特性作为参数输入，基于DEM 法对弯曲强度进行数值计算，再利用试验数据对其进行校核，验证数值计算方法的准确性，用以为模型冰力学特性预报提供参考。

1 基 于DEM 法 的 模 型 冰 弯 曲 强 度数值计算

1.1 DEM 法基本原理

1） 单元失效判断。

在离散元模型中，当内力平衡时，单元之间的相互作用可以看作是一个动态的过程，通过追踪单元的运动，可以得到每个单元的接触力和位移。考虑海冰的材料性质和离散特性，采用具有黏结−破碎特性的平行黏结单元对模型冰颗粒进行构造[7]。平行黏结模型是通过在两相邻颗粒接触区域添加附加材料，并设想有一组具有恒定法向和切向刚度的线性弹簧，以两个单元的接触点为中心，均匀分布在黏结圆盘上，从而完成单元间作用力和力矩的传递，如图1 所示。图中：Fn，Fs分别为接触力的法向力和切向力；Mn，Ms分别为法向力矩和切向力矩；xA，xB为单元中心点；R为黏结圆盘半径；n为法向方向[10]。

图1 单元之间的平行黏结模型Fig. 1 Parallel bonding model between elements

单元之间的黏结失效通过梁理论进行求解：

2） 冰层离散化。

在计算两个单元的接触力之前，需要确定单元之间的接触方式，判定其接触条件。本文首先对单元进行编号处理，采用元胞数组储存颗粒的位置和运动信息，并随着时间步的更新，利用元胞数组的特性对粒子的位置和运动信息进行实时更新，简化颗粒搜索过程，同时可实时计算粒子间的相互作用力。矩形单元采用最密排列的方式对冰层进行构造，冰层离散单元通过元胞数组进行编号（k1，k2，…，kn），然后利用冰单元kn计算相互间的作用力以及与其他冰单元之间的相互作用，如图2 所示。

图2 冰单元的编号方式Fig. 2 Numbering of ice units

然后，设置冰单元之间力的传递方式。根据冰层受力位置的不同，在冰单元之间力的传递也会有所不同。当左侧冰单元（k1，k2）受力时，冰单元之间的接触力会自左侧向右和向后传递；当右侧冰单元（k4，k5）受力时，接触力将自右侧向左和向后传递；当中部冰单元（k3）受力时，冰单元之间的接触力自中间向两侧和向后传递。完成第1 步的传递之后，接触力将逐渐拓展，直至覆盖到整个冰层。最后，通过对接触力的不同传递方式进行模拟，计算冰单元之间力的作用，以及冰层的裂纹生成与断裂。

3） 单元间受力分析。

在单元黏结处，单元受力和位移的关系可以通过以下参数来表示：法向刚度Kn和切向刚度Ks、摩擦系数µ、法向阻尼和切向阻尼，以及粒子重叠时在重叠区域中心形成的接触区。对接触处进行力学分析，接触力Fi可分解为法向力Fn和切向力Fs：

式中，ni和ti为单位向量。根据Hertz 接触理论，法向力Fn为弹簧和法向阻尼器作用在单元上的弹性力Fe=Knxn与阻尼力Fν=Cnvn的合力，可以表达为

式中：xn为单元之间的重叠量；Cn为单元间的法向阻尼系数；vn为法向相对速度。单元间的切向力Fs一般以增量的形式进行计算。剪应力增量∆Fs为每个时间步内切向方向的变化量与切向刚度的乘积：

式中：Δxs，vs分别为两个单元在接触点处的切向位移增量和相对切向速度；Cs为单元间的切向阻尼系数。假设单元间的接触力满足摩尔−库伦定律，则有

1.2 模型冰弯曲强度数值计算

在数值计算中，冰样的加载方式如图3 所示。图中，L为冰样长度，L0为加载点至冰样根部的距离，h为冰厚，b为冰样宽度，P为冰样上的作用力。在测试过程中，设置冰样根部为约束边界，其另一端为自由边界，以平行黏结圆盘黏结冰单元，以不同的加载速率在冰样端部向下施加载荷。冰单元的受力逐渐拓展至整个冰样，直至超过单元间的黏结强度，冰样发生弯曲破坏。根据悬臂梁的受力特点，冰样的弯曲强度 σf为

图3 用弯曲强度测试悬臂梁的方法示意图Fig. 3 Schematic diagram of a test method for cantilever beams using flexure strength

当冰样发生弯曲破坏时，弯曲测试力达到最大值Pmax，此时，计算得到的弯曲强度即为模型冰的弯曲强度 σf。

根据模型冰的力学参数试验数据，包括弯曲强度、压缩强度、弹性模量等，获得数值计算的细观参数如表1 所示。

表1 数值计算参数Table 1 Numerical calculation parameters

通过冰弯曲强度试验研究，确定数值计算工况，变换冰样的加载速率和冰样尺寸。其中，冰样加载速率的取值范围为1～5 mm/s，冰样宽度为2h～6h，冰样长度为4h～10h。此试验与冰的力学试验工况相同，旨在研究弯曲强度的加载速率相关性以及尺寸效应。确定了计算参数和工况之后，开展冰样弯曲强度的数值计算，图4 为加载速率为5 mm/s、冰样尺寸为240 mm×80 mm×40 mm时的计算时历曲线。

图4 弯曲强度计算时历曲线Fig. 4 Time history curve of flexure strength by calculation

2 模型冰弯曲强度试验测量

依托中国船舶科学研究中心小型冰水池（图5）的试验条件，针对盐水柱状模型冰开展弯曲强度试验[13]。该模型冰是在0.5% 的氯化钠溶液中通过自然冷冻形成。首先，将盐水混合物降温冷冻至冰点附近进行引晶，在水体表面形成一层极薄的细小颗粒冰晶层，进而在其下方引起以柱状结构为主的冰颗粒的生长；然后，在−18 ℃的条件下逐渐冷冻形成包含卤水泡的冰层，在达到目标厚度后，开始回温。在此过程中，冰层中的卤水慢慢流失，逐渐削弱冰层。达到目标强度后，维持温度不变，开始试验测试。

试验中，冰样厚度h的均值为40 mm。为研究冰样尺寸与弯曲强度之间的关系，冰样宽度取2h～6h，冰样长度取4h～10h。在测试过程中，以1～5 mm/s 的加载速率在冰样端部向下施加载荷直至冰样发生弯曲破坏，然后记录不同时刻冰样受到的作用力。试验时，每2～4 个冰样为一组，以保证结果的准确性。图6 所示为加载速率为5 mm/s、冰样尺寸为240 mm×80 mm×40 mm 时的加载过程和时历曲线。可见，冰样受力呈线性增加，冰样的破坏方式为脆性破坏。

3 数值计算与模型试验比较分析

冰的弯曲强度对冰区船舶破冰航行时的冰阻力有着重要影响，充分研究冰的弯曲性能，可以更好地了解冰的破坏方式及其与船体作用的原理。本文将通过一系列数值计算，确定冰样尺寸、加载速率和冰厚等因素的影响，并对数值计算模型进行校核，分析计算结果的误差率，以确保计算结果的可靠性。

将数值计算时历曲线与模型冰力学试验时历曲线进行对比，结果如图7 所示。由图可以发现，两者在冰样及弯曲过程中，时历曲线均表现为线性变化，力的增加较为均匀，冰样的破坏方式为明显的脆性破坏。

图7 数值计算与模型试验时历曲线对比Fig. 7 Comparison of time history curves between numerical calculation and model test

加载速率对冰弯曲强度的影响在不同的海域有着不同的表现。Timco 等[14]在统计分析大量的极地海冰研究数据后提出，加载速率对弯曲强度的影响并不明显。隋吉学等[15]通过在渤海的实测发现，加载速率与弯曲强度之间实际有着明显的相关性，且弯曲强度在不同的加载速率下有一定的韧脆性转变关系。而在统计模型冰试验数据后发现，弯曲强度与加载速率之间并无明显的相关性，弯曲强度趋于某一恒定值，如图8 所示。试验和计算结果与Timco 等[14]和Frederking 等[3]的海冰弯曲强度试验测量结果吻合，说明该模型冰与北极海冰的特性相似。对模型冰的破坏过程进行分析，显示在不同加载速率下，冰样晶体错位程度不同，冰样的破坏方式也不同。当加载速率较慢时，破坏裂纹逐渐生成，直至冰样破坏；当加载速率变快时，则表现为明显的脆性破坏。

图8 不同加载速率下的计算结果比较Fig. 8 Comparison of calculation results at different loading rates

如图9 所示，在试验中改变冰样长度后测量得到的模型冰弯曲强度其数据点十分离散，没有明显的规律性，冰样长度的改变与弯曲强度之间没有明显的相关性。在数值计算中，改变冰样的长厚比时，随着冰样长度的增加，发生弯曲破坏时冰样根部的转动幅度也会增加，导致冰晶颗粒错位充分，变形量较大；数值计算结果呈现出与弯曲强度试验结果相同的变化趋势，其间的相对误差均值约为12%。改变冰样的宽度之后，如图10所示，冰的弯曲强度随着冰样宽度的增加迅速增加，弯曲强度与冰样的宽厚比呈明显的线性相关，计算的相对误差均值约为15%，能够大致反映弯曲强度与冰样宽厚比之间的相关关系。

图9 冰样长厚比变化时计算结果比较Fig. 9 Comparison of calculation results when the ice sample's length-to-thickness ratio changes

图10 冰样宽厚比变化时计算结果比较Fig. 10 Comparison of calculation results when the ice sample's width-to-thickness ratio changes

通过对不同工况下数值计算结果与试验结果的比较发现，本文的数值计算方法能够模拟出弯曲强度随速度、冰样尺寸等参数变化的趋势，同时，也基本能够计算得到弯曲强度值，验证了计算模型的合理性。

4 结 论

冰的力学性质是研究冰区航行船舶冰阻力的重要参数。通过对模型冰力学特性的研究，可以为冰区航行船舶的设计、校核以及总冰阻力预报提供重要的参数依据。本文在利用原位悬臂梁方法对盐水柱状模型冰的弯曲强度进行测量的基础上，采用离散元模型对其进行了数值模拟。通过对计算结果的处理和分析，并与试验测量数据进行比较，得到以下主要结论：

1） 由于冰材料的特性较复杂，冰的弯曲强度受多种因素的影响，其中，冰样的宽厚比对其影响较大，弯曲强度基本不受加载速率和长厚比变化的影响。

2） 基于DEM 法建立的计算模型能够反映弯曲强度随参数变化而变化的趋势，得到了冰样弯曲强度值的基本范围，可为冰的力学特性预报提供参考。

在今后的研究中，将进一步考虑加载方式、加载方向以及浮力作用对弯曲强度的影响，从而系统研究模型冰的弯曲强度。