基于近场动力学的单晶冰弹性各向异性数值模拟方法1)

黄焱 王建平 孙剑桥

(天津大学水利工程仿真与安全国家重点实验室,天津 300350)

引言

近年来,环北极与近北极国家相继出台了新的极地政策或战略部署[1-3],极地海域的商业航运、资源开发、科学考察、环境保护及军事保障等活动也随之日趋白热化.海冰是极地海域标志性的环境条件,也是极地航行船舶和极地作业装备在安全性保障上面对的首要挑战.为应对这一挑战,学术界与工程界开展了大量关于冰-结构相互作用过程的数值模拟方法研究[4-6],以期能够为工程装备在现实冰环境条件下的承载过程提供可靠的预报.然而,天然海冰因其复杂的力学性质,致使当前的数值模拟结果直接面向工程应用时尚存在极大的阻碍.因此,如何建立可靠的冰力学性质的数值模拟方法,已成为当前学术界关注的焦点问题之一.

已有大量的研究表明,天然冰材料在变形和破坏行为上都具有显著的各向异性特征[7-10],这一特征成为冰-结构相互作用过程中产生复杂载荷过程的关键诱因.因此,在数值模拟中准确再现冰材料这一关键的力学特征,就成为保障数值预报结果合理性的必要条件.基于此,一些学者探索了将天然冰材料表现出的整体各向异性特征引入数值模拟方法的途径,如Castelnau 等[11]探索了一种模拟多晶冰黏塑性行为各向异性的自洽方法;Hibler 和Schulson 等[12]将缺陷导致的各向异性引入到海冰数值模型中,并据此进行了双轴压缩模拟,探讨了模型中海冰的裂纹扩展模式;季顺迎和岳前进[13]针对海冰的动力过程模拟问题,讨论了弹塑性各向异性模型的适用性;文献[14]基于原位冰芯样本的实测数据,建立了不同晶体织构下的弹性刚度张量,并将其引入冰层力学模型之中,由此计算和分析了多晶冰织构特征所导致的各向异性对冰内地震波传播方式的影响.但正如冰力学研究领域知名学者E.M.Schulson 所言,多晶冰材料(见图1)的各向异性特征根源于单晶冰的各向异性[15].所以,对于天然冰材料各向异性特征的模拟,首先应建立针对单晶冰各向异性的再现方法,由此构建出模拟天然冰材料各向异性特征的基本框架,再进一步由此框架扩展出多晶冰材料的模拟方法,才可能全面、准确地再现天然冰材料的特殊力学性质.本文将针对当前学术界缺乏单晶冰材料各向异性模拟方法的局面,开展探索性研究工作.

图1 S2 冰的水平向与垂向晶体切片[15]Fig.1 Horizontal and vertical sections of S2 ice[15]

另一方面,在数值模拟方法上,基于非局部作用思想的近场动力学方法,因其在统一连续介质力学理论与粒子模拟方法上具有突出的潜质,已在学术界引起广泛关注[16-19].在近场动力学模型中,物质点间的相互作用力函数(即“本构力函数”)包含了材料的物性信息,它不再以传统的应力-应变关系等形式出现,也不假设位移场的连续性,求解时不存在空间求导.此外,位移场的连续与否并不影响基本方程的求解,不连续性(如裂纹、界面)在求解运动方程时自然产生,因而不会出现传统连续性模型求解时出现的病态特征[18].当前,近场动力学方法已被大量应用于岩石、复合材料裂纹萌生扩展过程的模拟[20-22],并在材料复杂力学行为特征的模拟上显示出突出的优势.鉴于此,这一方法也开始被引入冰-结构相互作用问题的模拟中,并因其在冰材料断裂破坏现象模拟上的良好表现而备受青睐[23-26].需要指出的是,目前近场动力学方法在应用于工程预报时,都还是基于简单的各向同性材料假设,进而导致已有的模拟结果尚难以合理再现真实作用场景.

据此,本文在回顾单晶冰弹性各向异性力学表征的基础上,采用近场动力学方法建立单晶冰态基模型,通过对单轴压缩过程的数值模拟,探讨弹性各向异性特征的引入方式及相关的参数标定方法,以期建立合理可靠的单晶冰弹性各向异性数值模拟方法.

1 单晶冰弹性各向异性的力学特征

在小位移条件下,应力与应变之间的关系为

其中,σi为应力分量,i,j=1,2,···,6;εi为应变分量;Sij为柔度矩阵分量.

对于单晶冰所具备的六方晶系结构,如果假设其基面内符合力学各向同性,则可分别将晶体的C轴方向及其基面内的任意一对正交方向分别定义为X3,X1和X2应力主轴,相应的柔度矩阵中只有5 个非零参数[9]

Fletcher[9]收集整理了大量的单晶冰内弹性波传播的实验测试数据,并将这些数据回归为描述单晶冰弹性各向异性的柔度矩阵表达,矩阵中各项系数的值分别为:S11=1040 GPa−1;S12=−430 GPa−1;S13=−240 GPa−1;S33=850 GPa−1;S44=3140 GPa−1.

进一步利用张量变换可获得任一方向有效柔度的表达式为[9,15]

其中,Φ表示任一加载方向与C轴的夹角.而各个变形方向上的杨氏模量为该方向有效柔度的倒数

其中,Eref为参考杨氏模量.该值基于Fletcher 实验测试获得的单晶冰柔度矩阵系数计算得到,并将作为本文数值模型后续标定的参考值.

图2 为由式(3)和式(4)绘制得到的单晶冰弹性各向异性特征表征图:图中的实线描绘的是Eref在与C轴平行的任意平面内的空间分布情况,该曲线上任意一点的横纵坐标比值为该平面内变形方向与C轴夹角Φ的正切值,该点距离坐标原点的长度为相应方向上的杨氏模量值;图中的虚线表征以C轴方向杨氏模量为特征值的各向同性边界.

图2 单晶冰弹性各向异性表征图[9]Fig.2 Reference Young's modulus as a function of orientation[9]

由图2 可知,杨氏模量的分布是一条关于坐标原点中心对称和坐标轴双对称的闭合曲线,故对各向异性的描述可限定于第一象限内.

2 态基近场动力学模型

2.1 本构模型

在近场动力学模型的发展历经了两个阶段,即键基模型和态基模型.键基模型的理论基础会导致恒定泊松比的天然缺陷[18],为弥补这一理论缺陷,又进一步发展出了态基模型[27].因此,态基模型已成为当前近场动力学模型应用中的主角[28-29].下面对态基模型的基础理论进行简要介绍.

在态基近场动力学模型中,参考构型中x位置处材料点在t时刻的运动方程为

其中ρ为参考构型中材料的密度;u为位移向量场;Hx为x材料点的“近场范围”;T(x′,x,t) 为材料点x′施加在x上的力密度向量;T(x,x′,t) 为材料点x施加在x′上的力密度向量;T(x′,x,t)−T(x,x′,t) 对应于键基模型中的对点力函数;dVx′ 为材料点x′的体积;b为体积力向量场.

材料点x′和x之间的相互作用通过“键”来实现.以ξ表示参考构型中材料点x′和x的相对位置

在t时刻的当前构型中,材料点x′和x的相对位移为

故向量 ξ +η 表示两材料点在当前构型中的相对位置.

定义材料点x′和x之间“键”的伸长率为

其中,当“键”处于拉伸状态时,伸长率为正;当“键”处于压缩状态时,伸长率为负.

定义材料点x处的膨胀度为

其中d为与材料无关的参数.

在线性态基近场动力学固体模型中,力密度向量与当前构型中的相对位置向量 ξ +η 平行,表示为

其中,ω为影响函数,可根据材料点x′,x之间的距离控制它们的相互作用大小;a和b为与材料有关的参数;Λ 为

在近场动力学中,认为处于近场范围内的材料点均会与中心材料点发生相互作用,即非局部相互作用,这是近场动力学与传统连续介质力学最重要的区别之一.

文献[30]认为,影响函数应当包含 δ 作为长度尺度的对近场范围内其他材料点的影响,因此其提出以下影响函数

将影响函数代入式(10)可得态基近场动力学的力密度向量表达

其中,k和G分别为体积模量和剪切模量,在三维情况下分别为

其中,E为杨氏模量;ν 为泊松比.

2.2 时间积分格式

近场动力学运动方程是一个积分-微分方程,它的解通过空间和时间的数值积分来获得.本文采用无网格离散化方法,将材料离散为一系列具有确定体积的粒子,离散后粒子xi的运动方程变为

其中N为粒子xi近场范围内的粒子数,下标m表示第m个时间步.

粒子xi膨胀度θi为

其中,下标ij代表该物理量与粒子xi和xj相关.

在获得tm时刻的加速度后,采用以下显式积分算法获得下一时间步的速度和位移[30]

其中Δt为计算时间步长.

显式积分算法虽然简单明了,但是其只是条件稳定的,当时间步长大于最大容许时间步长之后便不能获得稳定收敛解.Madenci 等利用冯·诺依曼稳定性分析得到了临界时间步长为[30]

3 基于各向同性假设的单晶冰基础态基模型

本文的近场动力学模拟借助C++程序和Visual Studio 平台实现,有限元模拟采用ANSYS 软件完成.

在探究如何将各向异性引入单晶冰近场动力学模型的方法之前,首先以图2 中展示的各向同性假设边界为基准,对模型粒子数量(或粒子间距)取值进行标定.

如图2 所示,各向同性假设下的杨氏模量为11.76 GPa,相应地,体积模量k=8.7 GPa;剪切模量G=4.62 GPa.

不论是天然淡水湖冰还是海冰材料,其中包含的晶粒直径大多集中在2～ 10 mm 范围内[31-33].据此,本文采用的单晶冰试样模型尺寸为2.5 mm × 2.5 mm ×5.0 mm,密度为0.935 mg/mm3,如图3 所示.

图3 各向同性模型及载荷条件(单位:mm)Fig.3 Isotropic model and loading conditions (unit:mm)

令n为试件横截面长度与粒子间距的比值.本文中,n分别取为8,10,16,20,25,40,此时粒子间距 ∆x分别为0.3125 mm,0.25 mm,0.15625 mm,0.125 mm,0.1 mm,0.0625 mm,据此建立了6 种模型,相应的粒子数分别为1024,2000,8192,16000,31250,128000.近场范围半径δ取为3∆x[34-36].Silling 和Askari[34]建议将边界条件通过一个有限厚度带状区域内的粒子定义.本节向模型顶部和底部同时施加大小相同、方向相反的恒力压缩载荷,大小为F=100 N,通过厚度为δ的带状区域内的粒子体积力定义b

其中,M为顶部或者底部δ厚度区域内的粒子数.

由于本节施加的载荷为静力载荷,为获得稳定解,借助自适应动态松弛算法使模型快速达到稳定状态[37].同时采用商用有限元软件ANSYS 建立相同尺度和力学性质的有限元模型,施加相同载荷,并进行静力分析.提取近场动力学模型中位于Z轴上的所有粒子的Z向位移,与有限元计算结果进行对比,得到结果如图4 所示.

从图4 可以看出,当本模型中n从25 增加到40后,计算结果与有限元结果的差异均在2%以内,然而计算量却增加了近6 倍.因此,本文将 ∆x=L/25作为离散粒子间距,可以满足计算精度和效率需求.

图4 各粒子间距模型中心轴Z 向位移与有限元结果的对比Fig.4 Displacement variation along the loading direction

4 各向异性特征的引入方式与参数标定

在近场动力学理论中,要体现材料的弹性各向异性,需要从力密度的本构关系入手.在力密度向量的表达式(10)中,括号中含参数a的第一项表示体积膨胀度对力密度的贡献,含参数b的第二项可表示体积膨胀度以外的变形对力密度的贡献.本文认为,材料的弹性各向异性行为受体积变形和其他变形的综合影响,因此本文提出了基于影响函数ω的各向异性模型.

4.1 基于影响函数的弹性各向异性模型

Silling 等[27]在提出态基近场动力学理论时定义影响函数 ω 为:一个位于中心材料点近场范围内的严格为正的标量值函数.若影响函数仅与“键”长有关,即表示为,则影响函数为球形.目前,大量研究工作都是基于球形的影响函数[30,38-39],所表征的仅为材料的各向同性.Silling 等[27]指出,对于不规则形状边界的结构体或者混合物表面的模拟,可以通过定义形状不同的影响函数来实现.本文即遵循这一思路,拟通过定义合适的影响函数,以期使模型呈现出各向异性特征.

在一般的态基近场动力学理论中,力密度如式(10)所示.本文初步假定材料的各向异性完全由影响函数ωij控制,模型中的其他部分与各向同性材料相同,并且各向同性部分的材料参数以单晶冰C轴方向的力学属性为基准.

由于在各向同性理论中,力密度中的模量参数均与杨氏模量成正相关,因此,假设在各向异性模型中,各模量参数和杨氏模量具有相同的关于角度的表达,即力密度向量为

式(23)中的分母项可视为单晶冰杨氏模量随加载方向变化的归一化表征,由式(3)和式(4)得到;φ为模型中“键”与C轴的夹角.

模型约束端的总反力由所有与被约束粒子相关的力密度集体贡献,表示为

其中,R为总反力;M为被约束的粒子xi的数量;N为粒子xi近场范围内的粒子数.

在一般单轴压缩试验中,名义应力定义为两端反力与构件横截面积的比值,应变定义为构件标距段的伸长量与标距的比值,即

其中,F为反力;A0为构件变形前的横截面面积;∆l为标距段的伸长量;l0为标距.

类似地,定义近场动力学模型在单轴压缩模拟中的名义应力为底部反力R与变形前横截面面积的比值;标距为顶部虚粒子层底部与底部虚粒子层顶部之间的距离;名义应变为标距段的伸长量与标距的比值.结构如图5 所示.模型尺寸与第3 节中的基础模型相同,C轴始终位于XZ平面内.顶部的应变率载荷通过设置三层虚粒子向下的速度为应变率与试样高度的乘积来施加,底部固定约束通过刚性固定三层虚粒子的位置来施加.本文中应变率大小均为10 s−1.

图5 恒定应变率加载模型(单位:mm)Fig.5 Constant strain rate loading model (unit:mm)

在单晶冰中,C轴和垂直于C轴方向的力学属性是非常重要的两个属性;同时,单晶冰的杨氏模量在加载方向与C轴呈45°附近时最小,因此这三个方向的杨氏模量将是本文考虑的重点.在载荷沿Z方向加载的前提下,本文通过在XZ平面内顺时针旋转C轴的角度,实现C轴与加载方向夹角Φ的调节.采用式(23)的影响函数进行Φ=0°,45°,90°三个不同方向上的恒定应变率压缩变形模拟,得到名义应力随名义应变的变化关系.这三个加载方向的名义应力-名义应变关系如图6 所示.

图6 名义应力-名义应变曲线Fig.6 Nominal stress-nominal strain curve

对于线弹性模型,名义应力-名义应变曲线近似为直线,其斜率为杨氏模量.采用最小二乘法对该斜率进行拟合,即可得到数值模型的杨氏模量.在式(23)的影响函数下,三个加载方向的杨氏模量分别为:E0=8.461 GPa,E45=7.805 GPa,E90=7.877 GPa;下标数字表示加载方向与C轴方向的夹角.

观察模拟结果并将其与图2 的参考杨氏模量进行对比可以发现,尽管将杨氏模量与角度的关系直接代入影响函数中,能够使数值模型呈现一定程度的各向异性,但其结果与实验测试得到的参考杨氏模量存在较大差异,仍需对该影响函数进行修正.

4.2 影响函数修正方案一:引入辅助参数α

如前所述,影响函数虽然初步继承了杨氏模量随加载方向变化的力学表征关系,但计算结果与参考杨氏模量仍存在数量上的差异,尤其是沿C轴加载方向.据此,本文初步推断,若使模型能够模拟出单晶冰沿各个加载方向的杨氏模量,需要在影响函数中加入辅助参数α的控制

式中α可从整体水平上调节影响函数的值.式(23)可视为式(27)在α=1 时的特殊情况.此时的影响函数存在以下关系

式(28)表示当任意一根“键”ξij的影响函数中的α变为原来的ζ倍后,影响函数的值变为原来的1/ζ倍,在相同变形下,力密度也变为原来的1/ζ倍.根据式(25),反力也变为原来的1/ζ倍,则杨氏模量也应变为原来的1/ζ.

考虑到C轴方向的特殊性,首先调整α使沿C轴方向的杨氏模量与参考杨氏模量匹配.根据上一节计算得到的C轴方向的杨氏模量E0=8.461 GPa,令=0.7195.将该α控制下的影响函数代入数值模型,进行0°,45°和90°三个方向的加载模拟,最终得到该三个方向的杨氏模量分别为:E0=11.84 GPa,E45=10.81 GPa,E90=10.94 GPa.

将上述杨氏模量计算结果分别与α=1 时的计算结果进行对比,可以得到三个方向杨氏模量的比值分别为:0.7146,0.7220 和0.7200.这些比值与α=0.7195 非常接近,表明前述杨氏模量随α增大而减小的推断成立.

然而,尽管在α=0.7195 时沿C轴方向的数值模型杨氏模量与参考杨氏模量非常接近,但在另外两个方向的差异较大.进一步,绘制α=1 与α=0.7195下各加载方向的杨氏模量计算结果,并与参考杨氏模量进行对比,如图7 所示.可以发现,数值模型杨氏模量随加载方向变化的曲线形状在两种α值下几乎没有变化,即模型的各向异性程度没有得到改善.因此,仍需进一步修正影响函数中的角度控制参量.

图7 两种不同α 时的杨氏模量对比Fig.7 Comparison of Young's modulus for two different α

4.3 影响函数修正方案二:引入辅助参数α,β,γ

修正方案一的计算结果表明,仅通过单一的辅助控制参数并不能调整模型的各向异性程度.因此本文进一步推断,虽然影响函数与杨氏模量具有相似的关于角度的函数表达,但表达式中的各个分项的系数并不一定相同.据此,进一步假定影响函数形式为

其中,α,β,γ为独立变化的三个辅助控制参数.

将该影响函数代入力密度的公式,分别令φ=0°,45°和90°,得到沿这三个方向的力密度为

由式(30)可知,沿C轴方向“键”的力密度主要受α的控制,垂直于C轴方向主要受β的控制,其他方向将受到三个参数的综合影响.由于S44≥S11>S33,故与C轴呈45°夹角的“键”受γ影响最大.基于加载方向附近的“键”对该方向杨氏模量影响更大的思想,本文提出以下关于α,β,γ标定方法,上标(l)表示第l次标定

(1)l=0,令α(0),β(0),γ(0)的初始值均为0.7195,针对0°,45°,90°三个加载方向上的杨氏模量容许误差为Δ(本文设定为0.05%);

(2) 以α(l),β(l),γ(l)建立相应的近场动力学模型,进行三个方向的恒应变率单轴压缩数值模拟,分别得到对应的杨氏模量;

(3) 定义第l步时的总相对误差为

计算此时的总相对误差 ∆E(l);

(4) 若 ∆E(l)≤∆,输出α,β,γ的最终值为α(l),β(l),γ(l),反之,令

(5) 令l=l+1,重复第 (2)～ (4)步.

将α(0)=β(0)=γ(0)=0.7195 作为第一次循环的起始条件输入,执行上述标定过程.图8 展示了α,β,γ值以及三个方向上杨氏模量的相对误差随迭代次数的变化.可以看出,随着迭代次数的增加,α,β,γ分别稳定收敛于不同的值,同时三个方向的相对误差随着迭代次数的增加而下降,这表明本文提出的迭代标定方法基本有效.此外,总误差也随着迭代次数的增加迅速下降,如图9 所示.经过22 次循环后,总误差小于设定的容许误差Δ,此时α(22)=0.5364,β(22)=0.7506,γ(22)=1.2425,与此相对应的三个方向杨氏模量的计算结果分别为:E0=11.76 GPa,E45=8.793 GPa,E90=9.614 GPa.

图8 α,β,γ 及三个方向杨氏模量的相对误差随着迭代次数的变化Fig.8 α,β,γ and the relative error of Young's modulus as a function of iteration number

图9 总误差随着迭代次数的变化Fig.9 Total error as a function of iteration number

4.4 其他加载方向的适用性验证

影响函数修正方案二能够较为准确地模拟沿0°,45°和90°三个加载方向的弹性变形.进一步,对其他加载方向模拟结果的准确性开展验证.

将第22 次循环时的α,β,γ值作为终值,计算在此终值下与C轴方向分别呈15°,30°,60°和75°的四个加载方向的杨氏模量,及其与参考杨氏模量的相对误差,如表1 所示.图10 从曲线形状上展示了数值模型得到的杨氏模量与参考杨氏模量的对比.

图10 杨氏模量随加载方向的变化趋势对比Fig.10 Graphic comparison of Young's modulus under different loading directions

从表1 可以看出,其他四个加载方向的数值模型杨氏模量与参考杨氏模量的相对误差均位于5%以下,处于可接受的误差范围内.同时,数值模型杨氏模量随加载方向的变化趋势与参考杨氏模量较为契合.这进一步验证了,本文提出的基于影响函数的弹性各向异性近场动力学模型,能够实现单晶冰弹性各向异性力学表征的有效模拟.

表1 四个加载方向的杨氏模量计算结果及其相对误差Table 1 Young's modulus of the numerical model under different loading directions and the relative error

5 结论

本文基于态基近场动力学理论,根据单晶冰杨氏模量沿不同加载方向的变化规律,推断粒子间“键”的变形与力密度也存在类似的规律,提出了一种从影响函数出发的弹性各向异性模型.以前人实验测试得到的杨氏模量值为参考,通过开展与C轴呈0°,45°和90°三个加载方向的单晶冰单轴压缩数值模拟,提出了针对该影响函数的修正和辅助参数标定方法,并且最终在15°,30°,60°和75°等其他四个加载方向进行了验证.

计算结果表明,本文所提出的针对影响函数的修正与参数标定方法,能够较为便捷地找到数值模型杨氏模量与参考杨氏模量相一致的影响函数最优解,即本文提出的基于影响函数的近场动力学数值模拟方法,能够合理、准确地模拟单晶冰的弹性各向异性行为.

值得补充的是,极地海洋工程中常见的天然冰具有典型的多晶结构,其力学行为与晶体的排列方式、晶粒尺寸等有着密不可分的联系.同时,多晶冰的变形及裂纹成核与扩展模式不仅与单晶冰的变形有关,还受到晶界滑移等晶间作用机制的影响与控制.因此,在后续探索多晶冰的数值模拟方法时,仍需在单晶冰模拟的基础上,重点关注晶体织构的建模、单晶之间力的传递过程、断裂准则的引入等,以实现多晶冰变形与破坏各向异性的有效模拟.