高温冻土融化固结的水热力耦合模型构建1)

杨韬 郭颖 张程程 单炜

(东北林业大学，哈尔滨，150040)

随着全球气候变暖，多年冻土地区的平均地温升高，位于多年冻土南界及高海拔多年冻土外缘的不连续多年冻土分布区，普遍出现了退化现象[1-5]。温度的改变，使冰在冻土中的存在形式发生了变化，高温冻土中的未冻水含量的增加，使冻土由脆性逐渐向塑性转变[6-7]。因此，由于年平均地温的升高，冻土的稳定性变差，热敏感性提高，其上的建筑物及植被受到下部土体条件变化的影响也更为明显[8-9]。

研究人员最早结合太沙基(K. Terzaghi)一维固结理论，采用以0 ℃的温度界限为求解域边界的移动边界方法，求解冻土的一维融沉固结问题，但是在求解二维及三维融化固结问题时适用性较差[10]，Sykes et al.[11]通过结合三维固结理论及引入相变潜热作为温度场热源等方法，对其进行了改进。其后，Yao et al.[12]考虑到温度场参数在混合相变区是一个随温度变化的变量，并结合大应变理论，建立了三维大应变固结模型。Na et al.[13]应用有限应变下的广义硬化准则，结合线性运动平衡、物质平衡、能量平衡建立了水热力三项耦合模型，用以分析相变和传热作用下弹塑性响应和临界状态的演化过程。目前，对于冻土融化固结的研究，大部分采用的是以温度场求得的某一温度值为界限，温度场与流固耦合非同步的求解方式，只有少数研究考虑了冻土区变形的水热力三场同求解域求解[13]。

本研究依据传统的比奥(Biot)固结理论，以未冻水为耦合点，结合非稳态的对流传热方程，通过引入高温冻土应力-温度耦合损伤模型，建立水热力三个物理场同求解域的耦合模型，依据该模型分析高温冻土融化固结过程。

1 研究方法

1.1 耦合模型的建立

1.1.1 扩展的比奥大变形固结模型

比奥理论依据太沙基有效应力原理，建立了孔隙水压消散与土骨架变形的流固耦合关系[14]。冻土融化固结过程中，当施加荷载较大时，其即时构型与参考构型不再满足小应变状态时的近似重合假设[12]。因此，依据总拉格朗日(Total Lagrangian)方法建立大变形应力场控制方程，使模型可以更为准确地分析冻土融化固结过程中的几何非线性问题[15-16]。

在求解域内，张量形式的平衡方程可用式(1)表示。

∂Tij/∂Xj+F0i=0，i=1、2、3，j=1、2、3。

(1)

式中：Tij为第一类皮奥拉-基尔霍夫(Piola-Kirchhoff)应力；Xj为物质坐标；F0i为初始时刻参考构型之上的体积荷载。由于第一类皮奥拉-基尔霍夫应力张量为非对称张量，在后续计算中存在一定的不便。因此引入基矢量完全位于参考构型之上的第二类皮奥拉-基尔霍夫应力张量(Slk)对平衡关系进行分析，并将其转换为物质坐标的表达式。

(∂/∂Xk)[Slk(δil+∂ui/∂Xl)]+F0i=0。

(2)

(3)

引入不含空间坐标格林(Green)应变张量，建立宏观变形与微观应变之间的关系。

εij=(1/2)[(δki+∂uk/∂Xi)(δkj+∂uk/∂Xj)-δij]=

(1/2)[∂ui/∂Xj+∂uj/∂Xi+(∂uk/∂Xi)(∂uk/∂Xj)]。

(4)

εij表示的是宏观变形在微观应变上的总体映射，分析宏观变形发生的成因。εij由两部分组成，应力作用导致的土骨架变形(εeij)以及冰相变导致的体积收缩(εTij)，εTij可用式(5)表示。

εTij=(1/3)εVT=(1/3)(0.09/1.09)(dθice/dt)，i=j。

(5)

式中：εVT为冰相变引起的体积应变，作为一种等向变形，仅对应变球张量作出贡献；θice为体积含冰率，为该耦合模型的耦合节点。

通过广义胡克(Hooke)定律，建立应力应变两者之间的关系方程，即本构方程。

(6)

式中：Dijlm为本构张量。

依据质量守恒，孔隙内部的未冻水被排出，系统的质量发生了变化，则水力场的控制方程为：

-(k/yw)Δp=dεe/dt。

(7)

式中：k为渗透系数；Δ为二阶偏微分算子。

依据能量守恒，高温冻土单元温度场的控制方程如下：

Liρi(dθice/dt)+λpfΔT=Cpfρf(dT/dt)+

Cpwρw∑[∂(Tuj)/∂xj]。

(8)

式中：Li为冰的相变潜热；ρi为冰密度；λpf为导热系数；Cpf为高温冻土的恒压比热容；ρf为高温冻土的密度；Cpw为水的恒压比热容；uj为高温冻土形成的多孔介质内流体的流速。

1.1.2 耦合项

体积含冰率(θice)不仅对各物理场的因变量产生影响，而且成为了三个物理场共有的物理量，可由未冻水含量推导得出。

由于土颗粒表面能的存在，未冻水的含量与温度是一个非线性的关系，未冻水与温度的经验公式为[17]：

wu=α·|T|-β。

(9)

式中：α、β为试验系数；wu为未冻水含量。通过相关计算可得出wu与θice的关系：

θice=ρsρw(w0-wu)/[ρiρw+(ρi-ρw)wu+ρsρww0]。

(10)

式中：w0为初始含水率；联立式(9)、(10)即可得出θice与温度T的关系式，消除了控制方程的多余未知量。

1.2 模型参数计算方法

1.2.1 力学参数计算方法

依据损伤力学应变等价理论[17]，建立了高温冻土应力-温度耦合损伤本构模型，普通的杨氏模量被损伤模量所替代。

E′=E(1-Dc)。

(11)

式中：E为材料初始状态下的杨氏模量；Dc为复合损伤因子，可通过式(12)计算[18]：

Dc=DM+DT-DM·DT。

(12)

式中：DM为应力损伤因子；DT为温度损伤因子。

冻土受力劣化，应力损伤因子(DM)可通过某一时刻材料损伤单元数与初始总单元数的比值确定，因此，在某一应力水平下，损伤因子(DM)可表示为：

(13)

式中：Nt为损伤单元数；N为总单元数；η、F0是韦伯(Weibull)分布的参数。

依据莫尔-库伦(Mohr-Coulomb)强度准则，定义三维轴对称应力状态下的应力损伤表征量(F)[19]：

F=Eε1[σ1(1-sinφ)-σ3(1+sinφ)]/(σ1-2vσ3)。

(14)

式(13)中的2个试验参数(η、F0)，可通过直线拟合试验参数求解。由三轴压缩试验应力应变关系可得：

ln{-ln[(σ1-2vσ3)/Eε1]}=ηlnF-ηlnF0。

(15)

由上式通过拟合应力应变曲线，可直接得出参数η、F0的值。

随着温度升高，冻土内部冰相强度降低，胶结能力减弱。该种温度劣化效果，反映在冻土宏观层面即为冻土刚度降低，由此引入温度损伤因子(DT)进行分析。

DT=1-ET/ET0。

(16)

式中：ET为某一温度条件下对应的初始弹性模量；ET0为设无损状态下的初始弹性模量。弹性模量与温度的关系，可由式(17)表示[20]：

ET=β+γ|T|0.6。

(17)

联立式(11)、(12)、(13)、(16)，可得应力-温度耦合损伤模型。

1.2.2 水力及温度场参数计算方法

在水力场控制方程中，采用了表征土骨架透水能力的系数渗透率(k)。粉质黏土渗透率与孔隙比(e)的关系为：

k=CeD。

(18)

式中：C、D为材料参数[21]。

由于高温冻土中存在的冰相会对水分迁移形成阻滞作用，因此，冻结状态下高温冻土渗透率(kf)表达式为[22]：

kf=k/I=k/1010θi。

(19)

由于本研究的控制方程是依据高温冻土的多孔介质假设建立的，因此，使用复合材料的理论计算模型更为贴切，高温冻土的恒压比热(Cpf)为[12]：

(20)

式中：Cf为稳定冻结区的比热容；Cu为完全融化区的比热容；Tp、Tb分别为相变区的上下限温度(分别为1、-17 ℃)。

导热系数(λpf)为：

(21)

式中：λf为稳定冻结区的导热系数；λu为完全融化区导热系数。

2 结果与分析

2.1 验证试验方法

验证试验选取土样为哈尔滨市近郊粉质黏土(见图1)，相关物理力学参数：最大干密度1.86 g·cm-3、最佳含水率15.50%、液限32.70%、塑限19.50%。

图1 粒径级配累积曲线

采用该低液限粉质黏土，压实制备直径为10 cm、高为20 cm、含水率为21%的圆柱体试件，干密度为1.75 g/cm3。为减少试件冻结时产生的水分迁移量，将其置于低温试验箱内-25 ℃快速冻结，并恒定48 h。试验前24 h，置于恒温箱内-2 ℃，使试件整体温度分布均匀；随后置于低温三轴试验仪低温舱内，以-2 ℃恒温2 h，使系统温度恒定。

施加围压至0.25 MPa，设定目标轴向应力至1.57 kN；施加围压及轴压过程中，记录总轴向变形量。待达到目标轴压，且围压稳定后，以0.1 ℃/min速度均匀升高围压室温度；升温过程中，打开试件底部排水通道，即可实现对三向应力状态下高温冻结粉质黏土的融沉固结模拟试验。

2.2 模型参数计算

以-2 ℃为无损温度状态，通过0.05、0.15、0.25 MPa三组围压条件的低温冻土三轴压缩试验，确定应力损伤因子，试验曲线见图2。

图2 -2.0 ℃时三组围压条件对应的应力-应变曲线

依据力学参数计算方法对图2数据进行处理，可得应力损伤因子相关数据，围压0.25 MPa、E为141.64 MPa、F0为2.151 6、η为0.443 6、φ为14.66°。通过这些参数，即可得到围压为0.25 MPa时，应力-温度耦合损伤模型应力损伤因子。

在0.25 MPa围压时，-0.5、-1.0、-2.0 ℃三组温度条件对应的应力-应变曲线对模型温度损伤因子进行求解，试验曲线见图3。依据图3数据，拟合求解式(17)，最终结果见图4。融化状态下，则采用压缩试验所得压缩模量换算为杨氏模量予以确定，压缩模量(Es)为2.39 MPa。

图3 围压0.25 MPa时三组温度条件对应的应力-应变曲线

本研究粉质黏土融化阶段，存在一定的滞后性，其相变区上限为1 ℃[23]。因此，将式(9)进行变形，最终所得融化过程中未冻水含量与温度变化关系曲线见图5。

依据水力及温度场参数计算方法，确定水力及温度场参数。式(18)中C、D两参数，与材料的塑性指数、液性指数有关，表达式如下：

C=e-5.51-4ln(Ip)。

(22)

D=7.52e-0.25IL。

(23)式中：IP为塑性指数；IL为液性指数[21]。经计算，获得水力场、温度场相关参数(见表1)。

图4 围压0.25 MPa时弹性模量与温度关系曲线

图5 未冻水含量与温度关系曲线

表1 水力场与温度场相关参数

2.3 模型有效性检验

将计算获得的各项物理场参数带入水热力耦合(THM)模型控制方程进行计算，对比试验结果，验证模型有效性(见图6)。

由图6可见：A段内，高温冻土变形经历了缓慢变形及变形加速两个阶段，考虑几何非线性的THM模型计算结果，与试验曲线吻合度较高。B段为转变段，该处变形速率逐渐减小，模型计算结果与试验曲线出现少量偏离。C段为稳定段，此时模型计算结果与试验曲线在B段产生的偏离并未进一步扩大，两条变形曲线近乎平行变化。

A段为变形段；B段为转变段；C段为稳定段。

同时对比考虑几何非线性的大变形THM模型与不考虑几何非线性的THM模型，由图6可见：A段两者计算结果近乎没有差异；在B段，变形大于12 mm时，两者的差异出现。大变形模型的变形速率衰减速度明显小于小应变模型；但在C段，B段差异同样没有继续扩大，两者变形速率近乎相等。

融化固结为一个与时间相关的物理过程，为了分析变形速率的变化，定义试验固结度(U(t))，U(t)=uz(t)/uz(t1)，式中uz为一点竖直方向的位移、t1为最终试验结束时间，依此计算公式对沉降数据进行处理(见图7)。

图7 试验固结度试验曲线与模拟曲线对比

由图7可见：无论是考虑几何非线性的大变形THM模型，还是小应变假设下的THM模型，其固结进程实质上与试验曲线差别并不大。模型计算结果与试验结果固结曲线主要差别在于变形急速增加段后半部分，原因在于模型计算时，融土弹性模量采用了压缩模量换算而来的杨氏模量，其与冻结状态下的弹性模量存在不连续变化，造成了固结变形曲线的非平滑过渡。

2.4 融化固结过程分析

沿轴线方向取一垂直于上顶面与下底面的截面，选取0.9 h对比该截面上的温度及孔隙水压力(见图8)。该时刻试件温度处于高温冻土相变区域，可较好反映融化状态冻土内部物理场之间的相关关系。由图8可见：由于边缘先于中间部分融化，融化状态下的土体排水能力优于冻结状态，边缘土体形成排水通道，孔隙水压力先于中间部分消散；同时，孔隙水流动过程中，高温流体经过低温土骨架会交换部分热量，加速中间部分土体融化。

进一步对比应力场与温度场间的相关关系，同样选取0.9 h时的温度及范式等效(Von Mises)应力(见图9)。由图9可见：中部冻结部分土体所受应力高于边缘融化部分；由于冻土刚度高于融土，而上部施加轴向荷载的加压模块，迫使土体轴向发生整体变形，属于均匀变形，因此刚度较大的冻结部分承担了更多的应力，即为混合状态的土柱提供了更多刚度。

试验时间0.9 h；图面为孔隙水压力等值线图(云图)；图面数据为温度，单位为℃。

试验时间0.9 h；图面为范式等效应力等值线图(云图)；图面数据为温度，单位为℃。

最后对应力场及水力场间的相互关系进行分析，对比0.9 h时孔隙水压力和范式等效应力间的相关关系(见图10)。由图10可见：中部冻结部分应力最大，但是该处孔隙水压力并非最大；由于排水面位于冻土土柱下部，而位于顶端的区域孔隙水排出通道最长，因此孔隙水压力最大。依据有效应力原理，该处孔隙水压力承担了部分外荷载所造成的应力，则土骨架所受应力分布在该区域与孔隙水压力分布具有显著的相关性，随着孔隙水压力由边缘向中部逐渐增大，应力水平逐渐降低。

试验时间0.9 h；图面为范式等效应力等值线图(云图)；图面数据为孔隙水压力，单位为kPa。

3 结论

该模型通过引入高温冻土应力-温度耦合损伤本构模型，其冻结状态初始融化阶段预测结果与试验曲线较为吻合。在变形加速至缓慢变形过渡阶段(1.5～2.5 h)，模拟结果与试验曲线存在一定差异，但整体可有效预测高温冻结粉质黏土的融化固结过程。

考虑大应变状态的水热力耦合模型对融化变形的预测精度，优于依据小应变假设的模型。

高温冻土外部升温状态下，边缘融化部分形成排水通道，加速上部土体水分排出。