CN114578020A - 一种基于土体微观结构的正冻土体未冻水含量的计算方法 - Google Patents

Abstract

本发明涉及一种基于土体微观结构的正冻土体未冻水含量的计算方法，该方法包括以下步骤：(1)进行颗粒分析试验，得到试验土样颗粒的粒径分布曲线；(2)环刀取样后完全烘干，称得干土重量，继而得到土样干密度ρ_d；(3)得到该土样未冻水含量。本发明基于简化后的土体微观结构，结合毛管水弯液面半径和基质表面水膜厚度随温度的变化特征，推导出了毛管水体积含量的理论表达式，随后引入基于双电层理论计算结合水含量的简化表达式，将两处计算所得未冻水相加，得到总未冻水含量的理论公式。本发明需要的土样物性参数包括粒径分布曲线、干密度、比表面积，影响因子清晰，物理意义明确，可用于计算常规正冻土体任意温度下的体积未冻水含量。

Description

一种基于土体微观结构的正冻土体未冻水含量的计算方法

技术领域

本发明涉及冻土工程技术领域，尤其涉及基于土体微观结构的正冻土体未冻水含量计算方法。

背景技术

处于负温下的土体虽然早已达到自由水的相变温度，但由于土体基质的表面吸附作用及孔隙的毛细特性，仍保持有一定数量的液态水，称为未冻水。这些未冻水赋存在冰晶与土颗粒之间，对冻土的物理力学性质有很大的影响。一方面，冻结过程孔隙水持续相变成冰造成的体积膨胀效应会对周围土粒产生挤压作用，进而引起土体的冻胀现象，而气温回暖时孔隙冰的融化又会导致土体发生融沉现象，这两种现象极易对寒区工程建筑造成严重的损伤，甚至引起结构物的倒塌。另一方面，土体冻结时水分迁移会将盐分也带向冻结锋面，造成地表盐分的抬升，继而诱发地表次生盐渍化发育，无疑使脆弱的冻土区生态环境雪上加霜。因此，深入研究未冻水的存在机理，探索未冻水含量的影响因素对于指导寒区工程建设和加强生态环境治理具有重要意义。在研究土体的冻结融化过程中，进一步建立高精度的未冻水预测模型更有利于推动这一研究的发展，同时也能为数值计算奠定理论依据。

目前，对于正冻土体未冻水含量的计算方法使用较多的主要有以下两类：幂、指函数拟合类和土水曲线修正类。第一类中较为经典的为徐敩祖的经验模型，该模型形式简单，可清晰反应未冻水含量的变化，但缺乏必要的物理意义，需要用已有实测数据对模型未知参数进行拟合取值；第二类比较经典的有Zhen Liu提出的饱和度与温度的关系模型，这类模型一般是将“冻融过程”等同于“干湿过程”，继而引进土水特征曲线中较为成熟的表达式来描述冻结过程未冻水含量的变化，这类模型不仅形式复杂，且需额外拟合多个无物理意义的参数，局限性较大。

发明内容

本发明所要解决的技术问题是提供一种影响因子清晰、物理意义明确、计算精度较高的基于土体微观结构的正冻土体未冻水含量的计算方法。

为解决上述问题，本发明所述的基于土体微观结构的正冻土体未冻水含量的计算方法，包括以下步骤：

(1)进行颗粒分析试验，得到试验土样颗粒的粒径分布曲线；

(2)环刀取样后完全烘干，称得干土重量，继而得到土样干密度ρ_d；

(3)通过本发明计算方法，得到该土样未冻水含量：

负温下，土体孔隙中存在两种不同性质的水分，分别为毛管水和结合水。依据土粒与外侧水膜结合的牢固程度，又可将结合水细分为强结合水和弱结合水。现根据不同水分的性质特征，将土体内部孔隙结构进行简化；

根据土颗粒外侧水膜的分布情况，可取一个三维计算单元，投影后可得到平面计算单元。

平面计算单元中：r为毛管水弯液面的半径；2b为顶部和底部弯液面之间的距离；2a为两个相邻土颗粒最外侧结合水膜之间的距离；R为土颗粒的半径；h为两层结合水膜的总计厚度；毛管水顶部弯液面与最外侧结合水膜相交于P点，设P(x，y)；θ为毛管水弯液面与最外侧结合水膜的接触角，在土中，θ一般为0°；α为P点处弯液面的切线方向与水平方向之间的夹角，β为O′P与O′O之间的夹角。则各参数之间满足：

x²+(y-b-r)²＝r² (1)

[x-(R+h+a)]²+y²＝(R+h)² (2)

进一步地，可知：

x＝r·(R+h+α)·(R+h+r)^-1 (5)

式(5)中，只有当r＞a时得到的x值才可代入式(6)、(4)计算y和b；

为便于计算，可暂设P(c，d)，则图示平面计算单元中，自x＝0到x＝c的毛管水弯液面曲线绕x轴旋转一周所得旋转体体积V_tr可表示为：

化简后为：

此时可将式(8)中的c全部换为x，可得：

旋转体体积V_tr与土颗粒外侧结合水重合部分的体积V_coin可采用三重积分的形式表达，为简化计算过程，可取该部分体积独立建立坐标系；

故V_coin可表示为：

V_coin＝∫∫∫1dV (10)

积分区域Ω由XOY坐标平面与位于坐标平面上方的球面

组成。因而，上部分球面的z坐标可表示为：

将积分区域Ω投影至XOY坐标平面，则V_coin可表示为：

可进一步写为极坐标形式：

则：

此时可将式(14)中的d全部换为y，可得：

故三维计算单元中单个毛管水体积V_c为：

V_c＝V_tr-V_coin (16)

将式(9)、(15)代入式(16)，可得：

毛管水弯液面半径r(μm)与冰点降低值ΔT(℃)之间满足：

式中，σ为水的表面张力，N/m，其值与温度成线性关系，ρ_cw为毛管水密度，0.998g/cm³；T₀为273.15K；L是冰的融化潜热，336J/g；

孔隙率n、干密度ρ_d、土粒相对密度ρ_s之间满足：

孔隙率n可表示为：

结合式(19)、(20)，a值的大小可由下式确定：

结合水膜厚度h满足：

h＝δ+d (22)

弱结合水膜厚度d随温度T的变化关系可简化为：

式中，δ为强结合水膜厚度，一般可取

d为弱结合水膜厚度，单位为

T为负温的绝对值，单位为℃；

三维计算单元总体积为：

V＝8(R+h+a)³ (24)

三维计算单元中毛管水体积含量为：

故T_i温度时总毛管水体积含水率θ_c(T_i)可表示为：

式中，土颗粒半径R、该粒径所有颗粒体积占总体积的百分比v_d可由粒径分布得到；r可由式(18)得到；h可由式(22)、(23)得到；a可由式(21)得到；V_c可由式(17)得到。

T_i温度时总结合水体积含水率θ_b(T_i)可用比表面积A_s表示，并简化为：

θ_b(T_i)＝ρ_d·A_sh×10^-2 (27)

式中，A_s为土颗粒的比表面积，m²/g；h为结合水膜厚度，

ρ_d为土样的干密度，g/cm³。水的密度大小(g/cm³)接近于1，为便于计算，本式做了简化，所以并未体现出来，但等式右侧量纲为1。

因此，T_i温度时总体积未冻水含量θ_u可表示为：

θ_u＝θ_c(T_i)+θ_b(T_i) (28)

本发明基于简化后的土体微观结构，结合毛管水弯液面半径和基质表面水膜厚度随温度的变化特征，推导出了毛管水体积含量的理论表达式，随后引入了基于双电层理论的结合水含量计算表达式，将两处计算所得未冻水相加，得到总未冻水含量的理论公式。本发明需要的土样物性参数包括粒径分布曲线、干密度、比表面积，影响因子清晰，物理意义明确，可用于计算常规正冻土体任意温度下的体积未冻水含量。

本发明与现有技术相比具有以下优点：

1.本发明从毛管水与结合水的主要性质特征出发，简化了土体内部孔隙结构，首次考虑了基质外侧水膜厚度变化对毛管水含量的影响，并以三维计算单元中毛管水体积占比的形式建立了毛管水含量的理论表达式，随之引入基于双电层理论计算结合水含量的简化表达式，结合二者，最终给出了输入参数为粒径大小及其占比、干密度、比表面积的常规正冻土体未冻水含量的理论公式。

2.本发明可以计算任意常规正冻土体的未冻水含量，计算精度高，影响因子清晰，物理意义明确。

附图说明

下面结合附图对本发明的具体实施方式作进一步详细的说明。

图1为土体微观结构简化示意图。

图2为土颗粒外侧水膜示意图。

图3为毛管水简化示意图。

图4为三维计算单元。

图5为平面计算单元。

图6为V_coin计算过程。

图7为水-气交界面上的表面张力σ与温度T之间的关系。

图8为粉质粘土数据验证结果。

图9为莫玲粘土数据验证结果。

具体实施方式

基于土体微观结构的正冻土体未冻水含量的计算方法，包括以下步骤：

(1)进行颗粒分析试验，得到试验土样颗粒的粒径分布曲线；

(2)环刀取样后完全烘干，称得干土重量，继而得到土样干密度ρ_d；

(3)通过本发明计算方法，得到该土样未冻水含量。

负温下，土体孔隙中存在两种不同性质的水分，分别为毛管水和结合水。依据土粒与外侧水膜结合的牢固程度，又可将结合水分为强结合水和弱结合水。现根据不同水分的性质特征，将土体内部孔隙结构进行简化，如图1。土颗粒外侧水膜的分布情况如图2所示。毛管水的形状如图3所示。

从图1可取一个三维计算单元(如图4)，投影后可得到平面计算单元(如图5)。

平面计算单元中：r为毛管水弯液面的半径；2b为顶部和底部弯液面之间的距离；2a为两个相邻土颗粒最外侧结合水膜之间的距离；R为土颗粒的半径；h为两层结合水膜的总计厚度；毛管水顶部弯液面与最外侧结合水膜相交于P点，设P(x，y)；θ为毛管水弯液面与最外侧结合水膜的接触角，在土中，θ一般为0°；α为P点处弯液面的切线方向与水平方向之间的夹角，β为O′P与O′O之间的夹角。则各参数之间满足：

x²+(y-b-r)²＝r² (1)

[x-(R+h+a)]²+y²＝(R+h)² (2)

进一步地，可知：

x＝r·(R+h+a)·(R+h+r)^-1 (5)

式(5)中，只有当r＞a时得到的x值才可代入式(6)、(4)计算y和b。

为便于计算，可暂设P(c，d)，则图示平面计算单元中，自x＝0到x＝c的毛管水弯液面曲线绕x轴旋转一周所得旋转体体积V_tr可表示为：

化简后为：

此时可将式(8)中的c全部换为x，可得：

旋转体体积V_tr与土颗粒外侧结合水重合部分的体积V_coin可采用三重积分的形式表达，为简化计算过程，可取该部分体积独立建立坐标系，如图6所示

故V_coin可表示为：

V_coin＝∫∫∫1dV (10)

积分区域Ω由XOY坐标平面与位于坐标平面上方的球面

组成。因而，上部分球面的z坐标可表示为：

将积分区域Ω投影至XOY坐标平面，则V_coin可表示为：

可进一步写为极坐标形式：

则：

此时可将式(14)中的d全部换为y，可得：

故三维计算单元中单个毛管水体积V_c为：

V_c＝V_tr-V_coin (16)

将式(9)、(15)代入式(16)，可得：

毛管水弯液面半径r(μm)与冰点降低值ΔT(℃)之间满足：

式中，σ为水的表面张力，N/m，其值与温度成线性关系，如图7所示；ρ_cw为毛管水密度，0.998g/cm³；T₀为273.15K；L是冰的融化潜热，336J/g。

孔隙率n、干密度ρ_d、土粒相对密度ρ_s之间满足：

孔隙率n又可表示为：

结合式(19)、(20)，a值的大小可由下式确定：

结合水膜厚度h满足：

h＝δ+d (22)

弱结合水膜厚度d随温度T的变化关系可简化为：

式中，δ为强结合水膜厚度，一般可取

d为弱结合水膜厚度，

T为负温的绝对值，℃。

三维计算单元总体积为：

V＝8(R+h+a)³ (24)

三维计算单元中毛管水体积含量为：

故T_i温度时总毛管水体积含水率θ_c(T_i)可表示为：

式中，土颗粒半径R、该粒径所有颗粒体积占总体积的百分比v_d可由粒径分布得到；r可由式(18)得到；h可由式(22)、(23)得到；a可由式(21)得到；V_c可由式(17)得到。

T_i温度时总结合水体积含水率θ_b(T_i)可用比表面积A_s表示，并简化为：

θ_b(T_i)＝ρ_d·A_sh×10^-2 (27)

式中，A_s为土颗粒的比表面积，m²/g；h为结合水膜厚度，

ρ_d为土样的干密度，g/cm³。水的密度大小(g/cm³)接近于1，为便于计算本式做了简化，所以并未体现出来，但等式右侧量纲为1。

因此，T_i温度时总体积未冻水含量θ_u可表示为：

θ_u＝θ_c(T_i)+θ_b(T_i) (28)

模型结果验证

本发明选取两种不同质地的土壤进行验证：

(1)粉质粘土(取自青海北麓河，该土质的物性参数源于柴明堂的实测资料，粘粒、粉粒、砂粒均按美国农业部分类标准划分，如下表所示)：

塑限/％	液限/％	塑性指数	含盐量/％	粘粒/％	粉粒/％	砂粒/％
							18.6	36.7	17.94	0.0028	15.44	63.7	20.86

由图8可以看出，未冻水含量在整个负温区间内呈衰减趋势，当温度降到－6℃时，土体孔隙所持有的水分损失了近81％，至－16℃孔隙内约90％的水分已发生冻结，当温度降至－16℃以后，未冻水含量的变化逐渐趋于平缓。图8结果表明，由理论计算结果绘制而成的曲线能准确反映未冻水含量的变化过程，且理论计算结果与8个实测值数据的平均误差小于0.49％，说明该理论模型能高精度地实现对正冻土体未冻水含量的预测。

(2)莫玲粘土：(该土质的物性参数取自徐敩祖的实测资料，粘粒、粉粒、砂粒均按美国农业部分类标准划分，如下表所示)

如图9所示，分别给出了b、c、d、e四种不同初始含水率和干密度条件下(ρ_d依次等于1g/cm³、1.08g/cm³、1.26g/cm³、1.33g/cm³)莫玲粘土中未冻水含量随温度变化的实测值和理论计算曲线。b、c、d、e四种不同条件下莫玲粘土中未冻水含量的理论计算结果与实测值的平均误差分别为0.76％，0.82％，0.49％，0.44％，因而可知，对于不同初始含水率和干密度条件下的土体，本发明的计算结果仍然精确可靠。

综上，本发明所述计算方法适用性广，物理意义明确，计算精度高，可满足正冻土体未冻水含量的计算要求。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (3)

1.一种基于土体微观结构的正冻土体未冻水含量的计算方法，其特征在于，包括以下步骤：

(1)进行颗粒分析试验，得到试验土样颗粒的粒径分布曲线；

(2)环刀取样后完全烘干，称得干土重量，继而得到土样干密度ρ_d；

(3)通过本发明计算方法，得到该土样未冻水含量。

2.如权利要求1所述一种基于土体微观结构的正冻土体未冻水含量的计算方法，其特征在于：根据土颗粒外侧水膜的分布情况，可取一个三维计算单元，投影后可得到平面计算单元；

平面计算单元中：r为毛管水弯液面的半径；2b为顶部和底部弯液面之间的距离；2a为两个相邻土颗粒最外侧结合水膜之间的距离；R为土颗粒的半径；h为两层结合水膜的总计厚度；毛管水顶部弯液面与最外侧结合水膜相交于P点，设P(x，y)；θ为毛管水弯液面与最外侧结合水膜的接触角，在土中，θ一般为0°；α为P点处弯液面的切线方向与水平方向之间的夹角，β为O′P与O′0之间的夹角。则各参数之间满足：

x²+(y-b-r)²＝r² (1)

[x-(R+h+a)]²+y²＝(R+h)² (2)

进一步地，可知：

x＝r·(R+h+a)·(R+h+r)^-1 (5)

式(5)中，只有当r＞a时得到的x值才可代入式(6)、(4)计算y和b；

为便于计算，可暂设P(c，d)，则图示平面计算单元中，自x＝0到x＝c的毛管水弯液面曲线绕x轴旋转一周所得旋转体体积V_tr可表示为：

化简后为：

此时可将式(8)中的c全部换为x，可得：

旋转体体积V_tr与土颗粒外侧结合水重合部分的体积V_coin可采用三重积分的形式表达，为简化计算过程，可取该部分体积独立建立坐标系；

故V_coin可表示为：

V_coin＝∫∫∫1dV (10)

积分区域Ω由XOY坐标平面与位于坐标平面上方的球面

组成。因而，上部分球面的z坐标可表示为：

将积分区域Ω投影至XOY坐标平面，则V_coin可表示为：

可进一步写为极坐标形式：

则：

此时可将式(14)中的d全部换为y，可得：

故三维计算单元中单个毛管水体积V_c为：

V_c＝V_tr-V_coin (16)

将式(9)、(15)代入式(16)，可得：

毛管水弯液面半径r(μm)与冰点降低值ΔT(℃)之间满足：

式中，σ为水的表面张力，N/m，其值与温度成线性关系，ρ_cw为毛管水密度，0.998g/cm³；T₀为273.15K；L是冰的融化潜热，336J/g；

孔隙率n、干密度ρ_d、土粒相对密度ρ_s之间满足：

孔隙率n可表示为：

结合式(19)、(20)，a值的大小可由下式确定：

结合水膜厚度h满足：

h＝δ+d (22)

弱结合水膜厚度d随温度T的变化关系可简化为：

式中，δ为强结合水膜厚度，一般可取

d为弱结合水膜厚度，单位为

T为负温的绝对值，单位为℃；

三维计算单元总体积为：

V＝8(R+h+a)³ (24)

三维计算单元中毛管水体积含量为：

故温度为T_i时总毛管水体积含水率θ_c(T_i)可表示为：

式中，土颗粒半径R、该粒径所有颗粒体积占总体积的百分比v_d可由粒径分布得到；r可由式(18)得到；h可由式(22)、(23)得到；a可由式(21)得到；V_c可由式(17)得到。

温度为T_i时总结合水体积含水率θ_b(T_i)可用比表面积A_s表示，并简化为：

θ_b(T_i)＝ρ_d·A_sh×10^-2 (27)

式中，A_s为土颗粒的比表面积，m²/g；h为结合水膜厚度，

ρ_d为土样的干密度，g/cm³。水的密度大小(g/cm³)接近于1，为便于计算，本式做了简化，所以并未体现出来，但等式右侧量纲为1。

因此，T_i温度时总体积未冻水含量θ_u可表示为：

θ_u＝θ_c(T_i)+θ_b(T_i) (28)

3.根据权利要求2所述的一种基于土体微观结构的正冻土体未冻水含量的计算方法，其特征在于：负温下，土体孔隙中存在两种不同性质的水分，分别为毛管水和结合水，依据土粒与外侧水膜结合的牢固程度，又可将结合水细分为强结合水和弱结合水，根据不同水分的性质特征，将土体内部孔隙结构进行简化。

Images