寒区冰膨胀力的有限元分析

李 达，苏安双，孙颖娜，阎 烁

(1.黑龙江大学 水利电力学院，黑龙江 哈尔滨 150080；2．黑龙江省水利科学研究院，黑龙江 哈尔滨 150080)

在对寒区水工建筑物的影响研究中，多以研究冻土对水工建筑物冻害的研究为主，并且较为完善。魏锦等[1]总结了水工建筑物由于冻土的原因产生的冻害类别，以及其冻害机理，并简述了处理冻害的防护措施；蔡梅等[2]增加了在冻融循环的过程中，水工建筑物混凝土老化引起的冰害。而对于冰对水工建筑物的损害影响研究颇少，研究课题也大多围绕冰害防治方面；苏盛奎[3]就冰荷载问题进行了分析，从冰的力学性能、冰荷载的理论和计算公式以及计算实例上较为完善的阐述总结；张元有[4]就黑河市水工建筑物的冰害问题展开研究，提出利用聚苯乙烯硬质泡沫板作为保温材料减少静冰压力对水工建筑物的损害；徐秉衡[5]总结了以往对水工建筑物冰害问题的处理经验，并着重对工程设计和运行管理方面提出了有效治理冰害的措施。水工建筑物冰害的造成主要原因是冰膨胀力、动冰撞击力以及冰与水工建筑物界面间黏结力综合产生的。因此在水工建筑物的设计过程中常常依靠经验公式的计算和按照水工建筑物抗冰冻规范对冰压力的取值和计算进行设计，没有从分析冰害机理层面进行分析，以至于在复杂环境下的冰害问题上并没有得到广泛性的解决。

在信息技术高度发达的21世纪，利用计算机技术能够模拟越来越多真实的情况，解决实际问题，在冰膨胀数值分析的问题上，黄焱等[6]就冰膨胀力问题展开过数值分析，但是其主要内容主要针对就冰膨胀原理及约束条件对冰膨胀力的影响，其他因素并未详细阐述；刘荔铭[7]以ANSYS参数化语言对冰荷载进行数值分析，分析过程中，完善了对初始温度、温升率、冰盖厚度和边界条件四个方面进行单独分析，其反映温度场变化率的过程中，围绕的是冰盖厚度同温度的关系；另外李梦姌[8]借用ABAQUS有限元分析软件对冰荷载进行分析，除了对冰荷载影响因素单独分析外，还利用软件特点研究了具有裂纹的冰盖在水位下降过程中对水工建筑物产生的危害。在以上数值分析的过程中集中对冰荷载影响因素进行了重点分析，然而数值分析对实际工况的模拟还具有不完善性，以及借用有限元分析的手段来研究冰问题的精准性是需要论证的。本文在介绍冰膨胀力有限元分析的同时，着重进行数值分析。

1 冰温度膨胀力有限元分析计算

本文使用由ANSYS公司提供的ANSYS Academic Student(ANSYS学生版)以及图像，其中利用ANSYS WORKBENCH 19.2有限元软件对冰温度膨胀力(静冰压力)进行分析，ANSYS WORKBENCH有限元数值分析软件相较于ANSYS CLASSIC，具有更完善的图形软件接口，能够建立较为复杂的制图，而不出现丢线、丢面的状况；同时，在做网格划分、材料属性填充和耦合场计算等数值分析时，能够自动选取单元、优化划分、操作规范化和快捷化，有效的增强计算能力。

1.1 有限元模型仿真

根据气象资料，黑龙江省寒冷期极其长，多达5个月的封冻期，冬季室外温度多数在-20 ℃左右，加上黑龙江省河流众多，封冻期冰厚一般都在1 m左右，因此对于冰盖的模拟仿真，为确保具有普适性，建立长100 m、宽30 m、高1 m的冰模型，将模型划分网格，将模型沿厚度方向分10层,模型共划分27 840个单元,123 733个节点，如图1；另外添加冰材料属性，见表1。

图1 冰模型图和有限元分析模型图

弹性模量/MPa热传导系数/(W·m-1· ℃)比热/(J·kg-1· ℃)密度/(kg·m-3)热膨胀系数/(×10-6/ ℃)23000.612096.1+7.116 t(t为冰温)916.715.2

1.2 有限元模型仿真计算

温度场的变化和边界条件是导致冰膨胀力产生的直接因素，若要探究冰膨胀力的影响必须从这两方面共同入手。

1.2.1 稳态温度场计算

对稳态温度场计算是为了完成对冰膨胀力影响的最初始物质状态条件，其他因素的改变必须要以此为基础，尤其是初始温度改变时应重新计算稳态温度场，再加以分析。分析过程如下：

设置环境温度和上表面温度为-20 ℃，由于下表面与水直接接触设置为0 ℃，由此可以认为在模型内部形成了稳定的温度场，冰盖温度场分布，如图2，稳态场直线图，如图3，得知其分布状态是线性分布的，而且每一平面上的温度是相同的。

图2 冰盖稳态温度场分布效果图

图3 冰盖稳态温度场分布直线图

1.2.2 瞬态温度场计算

根据黑龙江省封冻期的河流情况，在3、4月份气温日较差最高，达到10 ℃～24 ℃，2 h内最高温差也可达到10 ℃，沿用上节稳态温度场结果，施加动态环境温度和上表面为在初始温度-20 ℃，2 h升高温度10 ℃，下表面与水直接接触设置为0 ℃的条件，由此可以认为在模型内部形成了动态的温度场，得到冰盖温度场的效果图，如图4所示，得知其分布状态是非线性分布的，同样每一平面上的温度是相同的，由图像可以看出，因温升作用对冰盖温度场的影响区域主要分布在冰盖上部15 cm左右的范围，另外冰盖内部温度场的温升率变化也是由外向里逐渐减弱，直到15 cm处变化率接近于0，其他厚度下的温度场与稳态温度场数值保持一致。

图4 2.0 h冰盖瞬态温度场分布效果图

1.2.3 瞬态应力场计算

在变化的温度场中，随着温度的上升，冰内部温度场的变化会导致冰体膨胀，由于边界条件的约束而产生膨胀力，所以对于稳态温度场，内部温度场恒定不变，所以不产生温度应力，只有瞬态温度场中才会产生温度应力。

直立墙约束是在冰盖和水工建筑物作用中最常见的约束，在数值建模中设置约束条件为假定冰盖与水工建筑物的接触一面做全约束，然后对整个冰盖长度和宽度方向做位移约束，厚度方向自由，在加载上面模型分析的初始温度-20 ℃、2 h升温10 ℃的温度场，得到应力场分布情况，如图5所示，分析得知，冰盖在厚度方向发生了较大的膨胀变形，并由于应力集中最大主应力发生在距离冰面1/3处，且在边角处达到最大值，最大主应力值为-0.271 67 MPa，为压应力。

图5 直立墙约束冰盖应力场分布效果图

2 结果分析

在以往的研究中，对有限单元法的适用性问题上的探究较少；就有限单元法利用有限元软件模拟计算的精准性，本文将以通过规范和经验公式的方式论证其可靠性。

(1)在北方地区对水工建筑物的设计和维护中，常常无法避免对冰工程的探索，根据大量的研究试验和工况监测，我国制定了中华人民共和国国家标准《水工建筑物抗冰冻设计规范》(GB/T 50662—2011)，用以对预防和设计寒区水工建筑物，本规范极具代表性和准确性，规范中对于静冰压力部分的数值规定如下，见表2，为了验证对静冰压力有限元分析准确性，应选取表中5种厚度的冰进行检验。

(2)徐秉衡[5]在对冰害的研究过程中，总结出最大静冰压力经验公式为：

Fd=137.3k·δ·Cδ·mt

(1)

式中：Fd为静冰压力，kN/m；137.3为应力常数，kPa;δ为冰厚，m；mt为与升温持续时间t有关的时间系数，一般天气，δ<50 cm，t=6 h，mt=1.0，δ>50 cm，2 d连续升温t=30 h，mt=1.82；k为综合影响系数，一般处取4.0～4.5，库面很大的平原水库取5.0，小型水库取3.5左右；Cδ为与冰厚有关的换算系数，见表3。

计算得静冰压力随冰厚变化，见表4。

表2 静冰压力Fd

注：1.对于库面狭小的水库，将表中静冰压力值乘以0.87后采用；对于库面开阔的大型平原水库，将表中静冰压力值乘以1.25后采用。2.冰厚取多年平均最大值。3.所列静冰压力值系水库在结冰期内水位基本不变情况下的压力，在此期间水位变动情况下的冰压力需做专门研究。4.静冰压力值按冰厚内插确定。

表3 换算系数Cδ

表4 经验公式计算静冰压力随冰厚变化表

(3)依次按(1)中对其他冰厚进行有限元分析，生成应力场分布效果图，如图6～图9，不同冰厚下最大主应力值、转换静冰压力值见表5。通过和规范数据、经验公式计算对比得出，见表6，有限元分析具有一定的适用性和准确性。

图6 0.4 m冰盖应力场分布效果图

图7 0.6 m冰盖应力场分布效果图

图8 0.8 m冰盖应力场分布效果图

图9 1.2 m冰盖应力场分布效果图

冰厚δ/m0.40.60.81.01.2最大主应力σ/MPa-0.208 80-0.232 38-0.247 96-0.271 67-0.289 37静冰压力Fd/(kN·m-1)-83.520-139.428-198.368-271.670-347.244

表6 规范、经验公式、数值分析法静冰压力随冰厚变化对比表

3 结 论

(1)本文通过总结冰荷载对水工建筑物破坏的原因以及破坏的各类研究，着重就冰温度膨胀力有限元数值分析适用性问题展开模拟试验和论证分析，分析得出在寒冷升温的气候下，冰盖会产生升温膨胀，由于约束作用，冰温度膨胀力的应力场随着瞬态温度场的变化而产生，数值分析得到1 m厚的冰盖在初始温度-20 ℃、2 h升温10 ℃的温度场，最大主应力值为-0.271 67 MPa，应力值较大，极可能造成水工建筑物的损伤和破坏。

(2)最后对不同冰厚下进行数值分析，并对其结果同国标规范和经验公式进行对比，事实证明有限元分析足够作为对日后性能劣化分析的工具，但是由于实际工况的复杂性，不能够仅借用有限元手段对工程设计的研究，仍需要利用规范设计—经验公式计算—简化试验—数值分析相结合的方法来解决问题。

(3)本文的数值分析及论证，一方面证明了有限元法数值分析可以充分作为工程研究的方法之一；另一方面，也提供了寒区水工建筑物在不同冰厚的情况下，冰对建筑的应力影响范围，可为日后水工建筑物的冰害研究、建造设计、防护维修提供有效的理论依据。