一种基于预熔理论的冻土未冻水含量确定方法

1.本发明涉及寒区冻土未冻水计算领域，尤其涉及一种基于预熔理论的冻土未冻水含量确定方法。


背景技术：

2.未冻水的含量影响孔隙冰的结构以及冻土中土壤颗粒和冰骨架的胶结强度，土体的内摩擦角随温度的降低而升高，负温下，土体的强度更大。水和冰在土壤中的体积因水冰相变而发生变化，土体中各组分的体积份额的变化也引起了土体导热系数、容积热容量的改变，这些相关变量对地质构造和岩土工程结构的热状态有显著影响，在预测寒区工程的热稳定性时，能有助于确定传热速率和冻融深度。未冻水含量控制着冻土的渗透性，在寒冷环境下土壤—水系统的水分输送过程中起着重要作用。在盐碱地区，未冻水携带着溶质由于低温吸力从未冻结区向冻结区迁移，随后在冻结缘的冷端形成的冰透镜体，同时由于地下水的补给使冰透镜体加厚，阻碍了水份迁移和溶质运输。
3.未冻水模型可以分为经验模型和理论模型。经验模型通过试验获取的数据并拟合曲线，模型中的参数物理意义不明确，且随工况的不同需要重新在实验室标定参数。理论模型中有具体数学推导过程，但模型表达式复杂难以应用于数值模拟计算，表达式中的参数更多且更难获取，综上所述，模型的表达式形式和参数数量及获取的难易程度限制了未冻水模型在数值模拟和寒区工程中的应用。


技术实现要素：

4.为克服上述现有技术存在的问题，研究一种涉及的参数较少，且已给出取值范围，表达式简单的模型及确定方法显得尤为重要。
5.本发明的目的在于提供一种模型表达式简单、计算精确的基于预熔理论的冻土未冻水含量确定方法，给出了未冻水含量的计算方式，并确定土壤中水膜厚度随掺杂水平的变化与温度的变化规律，只需要测量待测土壤的常用物理指标，如盐浓度，颗粒级配，孔隙比，即可计算得到土壤在不同温度下的未冻水含量。
6.本发明的目的是通过以下技术方案来实现的：一种基于预熔理论的冻土未冻水含量确定方法，包括以下步骤：
7.s1.建立水膜厚度与过冷温度之间的关系。
[0008][0009]
式中：
△
t表示过冷度，即初始冻结温度降低值，即与273.15k的偏差值；rg是气体常数；tm＝273.15k；n
im
是土颗粒表面的杂质浓度；ah为hamaker常数。
[0010]
s2.建立未冻水含量f
l
与水膜厚度关系的模型。
[0011][0012][0013]
式中，re表示颗粒的等效半径；d
p
表示冰颗粒和固体颗粒之间的水膜厚度；d
gb
表示冰颗粒之间的水膜厚度；sc和fcc表示两种常见堆积，简单立方堆积(simple cubic packing，简写为sc)和立方最密堆积(face-centered cubic packing，简写为fcc)，f
p
(sc)和f
p
(fcc)分别表示sc和fcc排列下颗粒占总体积的百分比(填充率)，可通过孔隙比e变换得来 (f
p
＝1/(1+e))；r表示冰液接触面的半径，通过gibbs-thomson关系确定。
[0014]
s3.采集土壤样本，测量土壤粒径分布，盐浓度，确定土壤中各粒径下土颗粒的体积百分比。
[0015]
s4.基于等效粒径表征土颗粒的粒径分布，对粉土、粉质粘土、黄土、砂土分别对等效粒径缩小0.21～0.23、0.16～0.28、0.25～0.3、0.36以修正土颗粒形状或表面粗糙度的影响；利用最大水膜厚度和已测得的盐浓度计算颗粒表面杂质密度。
[0016]
s5.计算冻土的未冻水含量。
[0017]
本发明所涉及的参数较少，且已给出取值范围，表达式简单。
附图说明
[0018]
下面结合附图对本发明的具体实施方式作进一步详细的说明。
[0019]
图1为本发明的方法流程图。
[0020]
图2为本发明基于等效粒径表征土颗粒粒径分布。
[0021]
图3为实施例sc模型；
[0022]
图4为实施例fcc模型中；
[0023]
图5a为水膜厚度与过冷温度的关系图之一；
[0024]
图5b为水膜厚度与过冷温度的关系图之二；
[0025]
图5c为水膜厚度与过冷温度的关系图之三；
[0026]
图6a为实施例第1种类型土壤的未冻水含量模型结果；
[0027]
图6b为实施例第2种类型土壤的未冻水含量模型结果；
[0028]
图6c为实施例第3种类型土壤的未冻水含量模型结果；
[0029]
图6d为实施例第4种类型土壤的未冻水含量模型结果；
[0030]
图6e为实施例第5种类型土壤的未冻水含量模型结果；
[0031]
图6f为实施例第6种类型土壤的未冻水含量模型结果；
[0032]
图6g为实施例第7种类型土壤的未冻水含量模型结果；
[0033]
图6h为实施例第8种类型土壤的未冻水含量模型结果。
具体实施方式
[0034]
结合实施例说明本发明的具体技术方案。如图1所示的流程，包括以下步骤：
[0035]
s1.建立水膜厚度与过冷温度之间的关系。
[0036][0037]
式中：
△
t表示过冷度，即初始冻结温度降低值(与273.15k的偏差值)。rg是气体常数； tm＝273.15k；n
im
是土颗粒表面的杂质浓度；孔隙溶液摩尔密度ρ
l
在计算中可认为是常数；qm为每摩尔的融化潜热。
[0038]
所述步骤s1包括：
[0039]
s101.预熔化是所有固体的普遍现象。当体系处于超低温时，体系中的间隙完全被冰晶填充。当温度上升后，与基质接触的冰的界面发生表面熔化，冰晶与冰晶之间发生晶界熔化。水分的存在形式有薄膜水和间隙水，此时依据介质形式可将系统划分为三层，层一是球状固体层、层二是宽度为d的预熔准液态层(水膜)，最后一层为冰晶层。若该体系内液体含有一定量杂质，系统的单位面积的gibbs自由能可以写成：
[0040]
g(t,p,ns,n
l
,n
im
)＝μs(t,p)ns+μ
l
(t,p)n
l
+μ
im
(t,p)n
im +rgt(nslnas+n
l
lna
l
+n
im
lna
im
)+g
interface
(d)
ꢀꢀꢀꢀ
(1)
[0041]
式中：t，p为温度和压力，μs,μ
l
,μ
im
分别为每摩尔固体、液体和杂质的化学势，μs,μ
l
,μ
im
分别为固体、液体和杂质的活度，在正常条件下，冰的活度为1，ns,n
l
,n
im
分别为固体、液体、杂质每单位面积的摩尔数，rg为气体常数，g
interface
为与水膜相关的界面自由能。固体(冰) 和液相(水)在凝固点处于平衡状态，其热力学平衡条件为：
[0042][0043]
在恒压下，可以确定每摩尔固体和液体的化学势差为：
[0044][0045]
式中：s
l
,ss分别为液体和固体的熵，在tm＝273.15k，相变潜热qm为常数，δt为t
m-t。g
interface
有两种贡献：
[0046][0047][0048]
式中：f
dis
(d)为色散力对界面自由能的贡献，f
elec
(d)为由于液层中离子捕获表面电荷的贡献，hamaker常数ah是比较范德华力相互作用强度的常用指标，qs为表面电荷密度，对于单一材料的固相和液相，范德华的相互作用是微弱的，当存在电作用力时，其占主导贡献。ε为液态水的相对介电常数，ε0真空介电常数，e为电子电荷量，kb为boltzmann常数，na为阿伏伽德罗常数，κ为常数，取7.237
×
107m-1/2
·
mol-1/2
。
[0049]
对于理想的稀释溶液，-lna
l
＝ρi/ρ
l
联立式(2)、式(3)和式(4)可得：
[0050]
[0051]
式中：qm为每摩尔的融化潜热，ρ
l
为液体的摩尔浓度，n
im
为单位面积的杂质数，由于冰中的杂质的偏析系数很小(约为10-6
)，因此认为预融液中杂质分布是均匀的。
[0052]
s102.为了得到冰—水界面曲率对冻结温度的影响，本实施例忽略冰表面能的各向异性，假设一个半径为r的球状冰晶完全被液态水包围。考虑到界面的面积不再是一个常数，将式 (1)改写为：
[0053]
g(t,p,ns,n
l
,n
im
)＝μs(t,p)ns+μ
l
(t,p)n
l
+μ
im
(t,p)n
im +rgt(n
s lnas+n
l lna
l
+n
im lna
im
)+4πr2γ
sl
[0054]
(7)
[0055]
对于曲率主半径为r1和r2的界面的情况，可以用式(2)计算热力学平衡时液体表面弯曲引起的温度下降。它们之间的关系如式(8)所示：
[0056][0057]
式中：冰晶半径r于ns相关(ns＝4πr3ρs/3)，ρs为冰的摩尔密度，γ
sl
为界面自由能。
[0058]
从图5a到图5c中可以看出，对于较大掺杂的体系，水膜厚度几乎完全由杂质效应影响，而在后续的计算中土壤是一种掺杂量较大的多孔介质，因此可对式(6)进行简化得式(9)：
[0059][0060]
s2.建立未冻水含量f
l
与水膜厚度关系的模型。
[0061][0062][0063]
式中，re表示颗粒的等效半径；d
p
表示冰颗粒和固体颗粒之间的水膜厚度；d
gb
表示冰颗粒之间的水膜厚度；f
p
是球体颗粒占总体积的比例，可通过孔隙比e变换得来；r表示冰液接触面的半径，可通过gibbs-thomson关系确定。为了简化吉布斯-汤姆逊方程，可以排除冰
‑ꢀ
水界面上较大的曲率半径。因此，式(8)表示为：
[0064][0065]
所述步骤s2包括：
[0066]
基于两种常见的堆积方式：简单立方堆积(simple cubic packing，简写为sc)和立方最密堆积(face-centeredcubic packing，简写为fcc)，建立含水量与颗粒表面和冰晶表面水膜、颗粒接触处裂隙和颗粒边界处的关系。
[0067]
(1)含水量与颗粒表面和冰晶表面水膜的关系
[0068]
图3b和图4b中黑色实线包围的截面为单元体sc和fcc的截面，取截面所对应的正方体为单元体。n
p,i
为球半径为ri的土壤单位体积的土壤颗粒数，对于sc模型和fcc模型n
p,i
分别为1/(8r
i3
)和√2/(8r
i3
)，f
p，i
为球体颗粒占总体积的比例与n
p,i
存在换算关系，即f
p，i
＝4πr
i3np,i
。则土颗粒表面积表示为：
[0069][0070]
对于sc模型，两个相邻冰晶之间的接触面积，等于正方形面积与土壤颗粒面积的差值。考虑到ri远小于ri，计算公式可以写成式(12a)：
[0071]sgb,i
(sc)≈(2ri+2ri)
2-π(ri+2ri)
2 (12a)
[0072]
对于fcc模型，每个颗粒边界的面积等于三角形面积与土壤颗粒面积的差。可由式(12b) 近似计算：
[0073][0074]
sc模型中每个土壤颗粒有3个晶界(图3a)，fcc模型中每个土壤颗粒有8个晶界(图 4a)。因此，单位体积内冰晶晶界贡献的整个表面积可以估算为：
[0075]sgb,i
(sc)＝3n
p,isgb,i
(sc) (13a)
[0076]sgb,i
(fcc)＝8n
p,isgb,i
(fcc) (13b)
[0077]
土壤颗粒和颗粒(冰)边界每单位体积的未冻水含量可以写成：
[0078][0079]
(2)曲率引起含水量的变化
[0080]
具有曲率的区域有两种，一种存在于两个相邻球体之间(图3b)，这是由由多个粒子组成的孔隙中形成冰晶造成的。另一个是冰晶边界和球体之间的边缘(图3c)，这是由相邻冰晶在颗粒表面上接触时形成的。
[0081]
对于sc模型，一个球体有六个相邻的球体，对于fcc模型，一个球体有12个相邻的球体。因此sc填料的每个球体有3个裂缝，fcc填料的每个球体有6个裂缝。当ri远小于ri时，每条裂缝的体积近似为2πr
iri2
，所以：
[0082][0083][0084]
如果忽略相邻球体和冰—水界面的曲率，则截面上的液态水(图3b)为：s
edge，i
＝2(r
i2-πr
i2
/4)。对于sc和fcc模型，单元体土颗粒对应的孔隙数(图3b和图4c)分别为1个和8个。sc 模型中，每个冰晶边界的总边长等于接触土壤颗粒的弧长(图3b)，4
×
2πri/4＝2πri。fcc 模型中每个颈部的总边等于三个接触土壤颗粒的弧长(图4c)，3
×
2πri/6＝πri。因此，可以用下列公式计算单位土壤体积颈部的总含水量。
[0085][0086]
[0087]
由半径为ri的土粒组成的液体总体积分数可归结为四种贡献的总和：
[0088][0089][0090]
考虑不同颗粒半径的土壤组成，可得土壤中的总液体分数如下:
[0091][0092]
式中：fi为粒径i土粒体积占总土粒体积的百分比，ri为粒径为i的土粒。采用等效粒径表征土颗粒粒径分布，将土体简化成等效粒径球状堆积体系。具体计算如公式(19)所示：
[0093][0094]
因此，结合式(18)和式(19)得式(20)，用于计算未冻水含量，依据s3中得讨论可对标定下式的适用孔隙比。
[0095][0096][0097]
s3.采集土壤样本，测量土壤粒径分布和盐浓度，确定土壤中各粒径下土颗粒的体积百分比。
[0098]
s4.基于等效粒径表征土颗粒的粒径分布，对粉土、粉质粘土、黄土、砂土分别对等效粒径缩小0.21～0.23、0.16～0.28、0.25～0.3、0.36以修正土颗粒形状或表面粗糙度的影响；利用最大水膜厚度和已测得的盐浓度计算颗粒表面杂质密度。
[0099]
所述步骤s4包括：
[0100]
s401.如图1所示，首先采用等效粒径表征土颗粒粒径分布，将土体简化成等效粒径球状堆积体系。具体计算如公式(19)所示。
[0101]
s402.对粉土、粉质粘土、黄土、砂土分别对等效粒径缩小0.21～0.23、0.16～0.28、0.25～0.3、0.36以修正土颗粒形状或表面粗糙度的影响。即s2中得α对粉土、粉质粘土、黄土、砂土分别取0.21～0.23、0.16～0.28、0.25～0.3、0.36。
[0102]
s403.天然状态下土壤中得颗粒并不一定相互接触，因此，无论哪种排列方式介质之间总存在间距2δ，堆积单位体系中有：
[0103][0104]
[0105]
土壤孔隙比表示土中孔隙的总体积与土颗粒总体积的比值，与f
p
存在换算关系：e＝1/f
p-1，土颗粒之间相互接触(即δ＝0)状态下利用式(21)可分别可求出e
sc
和e
fcc
的最小值为0.35和 0.9。
[0106][0107][0108]
随着温度无限接近冰点时，假定冰晶处于孔隙中心处，无限小可近似认为一点，此时冰晶具有最大水膜厚度。
[0109][0110][0111]
依据颗粒排列几何关系建立不同土颗粒堆积下的孔隙初始体积浓度与水膜表面杂质浓度之间的关联，如公式(24)所示。
[0112]nim
＝ηc0d0ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(24)
[0113]
式中：η为单位分子杂质电解离子数，例如：氯化钠η＝2，c0为溶液的浓度。
[0114]
s5.计算冻土的未冻水含量。
[0115]
所述步骤s5可参考图1，具体步骤如下：
[0116]
s501.结合s3中测得的土壤参数并结合s4计算结果代入s1中已建立的水膜厚度与过冷温度关系的模型中。
[0117]
s502.根据土壤的实测孔隙比选取适合的堆积模型并结合水膜厚度与过冷温度的关系建立冻土的未冻水含量模型。
[0118]
图6a到图6h为8种类型土壤的未冻水含量模型结果。