冻结过程中单孔冰结晶变形机制研究

doi:10.7522/j.issn.1000-0240.2022.0129

冻结过程中单孔冰结晶变形机制研究

张永安^,¹, 应赛^,²^,³, 文桃²^,³, 王月礼⁴, 周自强³, 卢千禧², 张支立²

1.甘肃有色工程勘察设计研究有限公司，甘肃兰州 730000

2.长江师范学院建筑物全生命周期健康检测与灾害防治工程研究中心，重庆 408100

3.甘肃省科学院地质自然灾害防治研究所，甘肃兰州 730000

4.山西冶金岩土工程勘察有限公司，山西太原 030000

The study of deformation mechanism of ice crystals in a freezing single pore

ZHANG Yong’an^,¹, YING Sai^,²^,³, WEN Tao²^,³, WANG Yueli⁴, ZHOU Ziqiang³, LU Qianxi², ZHANG Zhili²

1.Lanzhou Non-ferrous Engineering Investigation Design Limited Company，Lanzhou 730000，China

2.Engineering Research Center for Health Monitoring in Building Life Cycle and Disaster Prevention，Yangtze Normal University，Chongqing 408100，China

3.Gansu Academy of Sciences，Lanzhou 730000，China

4.Shanxi Metallurgical Rock-Soil Engineering Investigation Limited Company，Taiyuan 030000，China

通讯作者: 应赛，讲师，主要从事寒区岩土工程和特殊土力学研究. E-mail： yingsai35910@outlook.com

收稿日期: 2022-05-01 修回日期: 2022-09-27

基金资助:

第二次青藏高原综合科学考察研究项目.  2019QZKK0902
重庆市教委科技项目.  KJQN202101443
甘肃省公益性行业科研专项资金子项目.  2017121002
甘肃省自然科学基金项目.  20JR10RA472
重庆市自然科学基金项目.  cstc2019jcyj-msxmX0813
重庆市教委科技项目.  KJQN202101447

Received: 2022-05-01 Revised: 2022-09-27

作者简介 About authors

张永安，高级工程师，主要从事水文地质、工程地质与环境地质勘查与灾害防治研究.E-mail：zha0301@hotmail.com , E-mail：zha0301@hotmail.com

摘要

孔隙水在冻结过程中产生的冰结晶压力导致了多孔材料的冻胀及破坏。本文通过理论分析分别给出了不同形状晶体的结晶压力计算模型，并分析了经典结晶压力计算公式的使用条件。建立了降温过程中孔隙冰晶生长模型，实现冰晶生长过程中的孔隙变形计算，分析了晶核密度、孔径大小、荷载和冰晶体积对孔隙冻胀变形的影响机制。结果表明：起始孔隙直径和长宽比的增大对结晶变形抑制作用的机理在于减少了冰晶体积中膨胀结晶的比例。荷载对孔隙变形的抑制机制在于，荷载的增大迫使冰晶更多地横向生长（长宽比增大），导致膨胀结晶所占比例减小。孔隙中的晶体生长有完全填充模式和部分填充模式，在部分填充模式下，晶核密度、荷载和孔径的增大都会导致晶体在孔隙中的填充率增大，从而对孔隙结晶变形产生影响。本模型揭示了单个孔隙中冻胀变形机制，为解决多孔介质的冻胀变形与破坏问题提供了新的思路。

关键词： 结晶压力 ; 变形机制 ; 冰晶生长

Abstract

The crystallization pressure on the surface of the pore wall lead to the frost heave of porous materials during the growth of ice crystals in the cooling process. In this paper， the calculation models of crystallization pressure for crystals of different shapes are given respectively， and the application conditions of the classical calculation formula of crystallization pressure are analyzed. Then， based on crystallization kinetics theory and considering the interaction between crystal and pore wall， a growth model of pore ice crystallization during the cooling process was established to calculate the pore deformation during the growth of ice crystals， and the influence mechanism of crystal nucleus density， pore size， load and freezing amount on the pore crystallization deformation was analyzed. The results show that the increase of initial pore diameter and aspect ratio inhibits crystallization deformation by decreasing the proportion of expanded crystals in ice crystal volume. The inhibition mechanism of load on pore deformation is that the increase of load forces salt crystals to grow more laterally （the length-to-width ratio increases）， leading to the decrease of the proportion of expansive crystals.There are two modes of crystal growth in the pore： full filling mode and partial filling mode. In the partial filling mode， the increase of crystal nucleus density， load and pore size will lead to the increase of crystal filling rate in the pore， thus affecting the pore crystallization deformation. This model reveals the deformation mechanism of ice crystals in a freezing single pore， which provides a new idea for solving the problem of frost heave deformation and failure of porous media.

Keywords： crystallization pressure ; deformation mechanism ; growth of ice crystals

PDF (2527KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

张永安, 应赛, 文桃, 王月礼, 周自强, 卢千禧, 张支立. 冻结过程中单孔冰结晶变形机制研究[J]. 冰川冻土, 2022, 44(5): 1429-1439 doi:10.7522/j.issn.1000-0240.2022.0129

ZHANG Yong’an, YING Sai, WEN Tao, WANG Yueli, ZHOU Ziqiang, LU Qianxi, ZHANG Zhili. The study of deformation mechanism of ice crystals in a freezing single pore[J]. Journal of Glaciology and Geocryology, 2022, 44(5): 1429-1439 doi:10.7522/j.issn.1000-0240.2022.0129

0 引言

降温过程中冰晶生长在孔壁表面产生的结晶压力是导致多孔材料产生冻胀变形及破坏的主要原因。冰晶与孔壁相互作用的微观机制，是学者们关注的焦点问题。针对孔隙水的冻结相变，赖远明等^［1-2］提出了温度驱动的孔隙水冻结与融化过程的相变速率方程。针对结晶压力，Scherer^［3］在微观层面上建立的孔隙结晶热力学模型，主要考虑了溶液浓度和晶体表面曲率对结晶压力的影响，并在岩石、混凝土、黏土砖及盐渍土等多孔材料结晶变形及破坏的研究中广泛运用^［4-8］。许多学者认为土中水结冰体积要膨胀8.8%，是原位冻胀的产生机理。而上述对冻胀机制的描述过于简单，只强调了在冻结发生时液体水向固态冰转化过程中产生了体积膨胀，无法清晰描述冻胀过程。

本文基于结晶动力学理论，考虑了孔隙中冰晶生长过程，建立了降温过程中孔隙冰晶生长模型，实现冰晶生长过程中的单孔变形计算，分析了晶核密度、孔径大小、荷载和结冰量对单孔冻胀变形的影响机制，为解决孔隙中的冻胀变形与破坏问题提供了新的思路。

1 冰结晶压力

1.1　大孔中的冰结晶压力

考虑冰晶直径大于1 μm的情况，此时冰晶存在于大孔中，可忽略冰晶曲率的影响，以标准状态为基准点（参考点），溶液中水和冰的化学势分别为：

μ_{l} (T, p_{l}, a_{w}) = μ_{l}^{0} (T_{0}, p_{0}) + \int_{p_{0}}^{p_{l}} V_{w} d p - \int_{T_{0}}^{T} S_{w} d T + R T l n a_{w}

（1）

μ_{i} (T, p_{i}) = μ_{i}^{0} (T_{0}, p_{0}) + \int_{p_{0}}^{p_{i}} V_{i} d p - \int_{T_{0}}^{T} S_{i} d T

（2）

式中：μ_i和μ_l分别为冰和溶液中水的化学势；p_i和p_l分别为冰和溶液的压力（Pa）；a_w为溶液中水的活度；T为温度（K）； $μ_{l}^{0} (T_{0}, p_{0})$ 和 $μ_{i}^{0} (T_{0}, p_{0})$ 分别为溶液中水和冰在标准状态的化学势（T₀=273.15 K，p₀为标准大气压）；V_w和V_i分别为水和冰的摩尔体积（m³·mol^-1）；S_w和S_i分别为水和冰的摩尔焓（J·mol^-1·K^-1）；R为理想气体常数（J·mol^-1·K^-1）。

冰-水共存体系中当冰与溶液中的水达到平衡时，二者的化学势相等，有

\int_{p_{0}}^{p_{l}} V_{w} d p - \int_{p_{0}}^{p_{i}} V_{i} d p + R T l n a_{w} = \int_{T_{0}}^{T} S_{w} - S_{i} d T

（3）

当标准状态下，孔隙溶液压力与大气压力平衡，有p₀=p_l，考虑晶体与孔壁接触，产生的结晶压力 $Δ p = p_{i} - p_{l} > 0$ ，代入公式（3）可得

Δ p = \frac{Δ S_{m}}{V_{i}} (T_{0} - T) + \frac{R T}{V_{i}} l n a_{w}

（4）

式中：△S_m=S_w-S_i=L_wi/T₀为冰水相变时的焓变（J·mol^-1·K^-1），L_wi=6 010 J·mol^-1为冰水相变潜热。

由于溶液中水的活度a_w与溶液的冻结温度T_f相关，公式（4）隐含了冻结温度T_f对结晶压力的影响，自由溶液的冻结温度可由下式计算^［9］：

T_{f} - T_{0} = \frac{R T l n a_{w}}{Δ S_{m}}

（5）

将冻结温度T_f代入公式（4）可得以下形式：

Δ p = \frac{Δ S_{m}}{V_{i}} (T_{0} - T_{f}) + \frac{Δ S_{m}}{V_{i}} (T_{f} - T) + \frac{R T}{V_{i}} l n a_{w}

（6）

将公式 $T_{f} - T_{0} = \frac{R T_{f} l n a_{w}}{Δ S_{m}}$ 代入上式可得：

Δ p = \frac{R T}{V_{i}} l n a_{w} - \frac{R T_{f}}{V_{i}} l n a_{w f} + \frac{Δ S_{m}}{V_{i}} (T_{f} - T)

（7）

式中：a_wf为T_f温度下与冰晶平衡时的溶液中水的活度。考虑温度由T_f变化至T，该过程中溶液浓度不变或变化十分微小的情况，此时有a_w≈a_wf，代入公式（7）可得：

Δ p = \frac{R (T - T_{f})}{V_{i}} l n a_{w} + \frac{Δ S_{m}}{V_{i}} (T_{f} - T)

（8）

从公式（8）可以看出，冰结晶压力不仅与温度有关，还与溶液的活度有关。当溶液为纯水时，有 $l n a_{w} = 0$ ，代入公式（8），得到纯水状态下的冰晶压力，与Saetersdal^［10］和Kurylyk等^［11］提出的一类冻结土Clapeyron方程的冰晶压力一致。

2.2　小孔中的冰结晶压力

当孔隙直径小于1 μm时，孔隙中冰晶的化学势受其表面曲率的影响，晶体表面曲率的增大导致晶体化学势和溶解度的增大^［12］，此时冰的化学势为：

μ_{i} (T, p_{i}) = μ_{i}^{0} (T_{0}, p_{0}) + γ_{s l} V_{i} \frac{d A}{d V} + \int_{p_{0}}^{p_{i}} V_{i} d p - \int_{T_{0}}^{T} S_{i} d T

（9）

式中： $γ_{s l}$ 为冰与溶液的界面能，A为冰晶表面积，V为冰晶体积。

利用公式（9）代替公式（2），经1.1节的推导过程可得，考虑交界面曲率的冰晶压力计算公式：

Δ p = \frac{R (T - T_{f})}{V_{i}} l n a_{w} + \frac{Δ S_{m}}{V_{i}} (T_{f} - T) - γ_{s l} \frac{d A}{d V}

（10）

上式说明冰晶压力受溶液中水的活度、温度和冰晶表面曲率的共同影响。当晶体的不同面上的曲率存在差异时，不同面上的结晶压力是不同的，利用公式（10），代入某个面上的曲率值，可得该面上的结晶压力。

当考虑孔隙表面曲率的影响时，孔隙中溶液的冻结温度 $T_{f}^{*}$ 要小于自由溶液冻结温度 $T_{f}$ ，公式（10）可以写成如下形式：

\begin{array}{l} Δ p = \frac{R (T - T_{f}^{*})}{V_{i}} l n a_{w} - \frac{R (T_{f} - T_{f}^{*})}{V_{i}} l n a_{w} \\ + \frac{Δ S_{m}}{V_{i}} (T_{f} - T_{f}^{*}) + \frac{Δ S_{m}}{V_{i}} (T_{f}^{*} - T_{f}) - γ_{s l} \frac{d A}{d V} \end{array}

（11）

文献［9］给出了 ${T_{f}}^{*}$ 与 $T_{f}$ 之间的关系：

T_{f} - {T_{f}}^{*} = \frac{R (T_{f} - {T_{f}}^{*}) l n a_{w}}{Δ S_{m}} + \frac{γ_{s l} V_{i}}{Δ S_{m}} \frac{d A}{d V}

（12）

将上式代入公式（11）可得以孔隙溶液中冻结温度 ${T_{f}}^{*}$ 为参量的冰晶压力表达式：

Δ p = \frac{R (T - {T_{f}}^{*})}{V_{i}} l n a_{w} + \frac{Δ S_{m}}{V_{i}} ({T_{f}}^{*} - T)

（13）

从上述可以看出，当 $T \geq {T_{f}}^{*}$ 时， $Δ p \leq 0$ ，此时冰晶压力为0，只有当 $T < {T_{f}}^{*}$ 时， $Δ p > 0$ ，冰晶压力开始产生。上式表明冰晶压力受溶液中水的活度、温度和冰晶表面曲率的共同影响，其中晶体表面曲率的影响，体现在 ${T_{f}}^{*}$ 这一项上。当孔隙直径大于或等于1 μm时，孔隙中冰晶的化学势受其表面曲率的影响很小，可以忽略，孔隙中的结晶压力可由公式（8）计算。需要着重指出的是，上述公式的推导是建立在冰水两相化学势平衡的基础之上，所以上述公式的应用范围仅限于存在冰水两相平衡的系统，若孔隙中的水完全冻结成冰时，此时冰水两相平衡不存在，上述公式不再适用。

2.3　球形大孔与小孔连接时的冰结晶压力

Everett^［13］给出的冰晶压力计算的经典公式，其中考虑了晶体表面曲率的影响，形式上与公式（10）有差异，但本质上是公式（10）的一个特殊形式，即球形大孔与小孔连接时的球形大孔中产生的冰晶压力，以下给出详细分析。

考虑一个球形大孔（孔径r₁）和一个小孔（孔径r₂）连接的情况，系统温度T等于小孔的冻结温度 ${T_{f 2}}^{*}$ ，大孔完全被冰晶填充，小孔中处于冰水平衡状态，如图1所示。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 球形大孔中晶体的生长示意图

Fig.1 Schematic diagram of crystal growth in a spherical pore

由于大孔中的冻结温度 ${T_{f 1}}^{*}$ 高于小孔中的冻结温度 ${T_{f 2}}^{*}$ ，此时大孔处于过冷状态，当小孔中的水进入大孔，供给大孔冰晶的生长，在大孔中产生的冰晶压力为（孔隙溶液为纯水）：

Δ p = \frac{Δ S_{m}}{V_{i}} (T_{0} - {T_{f 2}}^{*}) - γ_{s l} \frac{d A_{1}}{d V_{1}}

（14）

式中：A₁为大孔中冰晶表面积；V₁为大孔中冰晶体积。

由公式（12）可得：

T_{0} - {T_{f 2}}^{*} = \frac{γ_{s l} V_{i}}{Δ S_{m}} \frac{d A_{2}}{d V_{2}}

（15）

式中：A₂为小孔中冰晶表面积；V₂为小孔中冰晶体积。将公式（15）代入公式（14），可得：

Δ p = γ_{s l} (\frac{d A_{2}}{d V_{2}} - \frac{d A_{1}}{d V_{1}})

（16）

考虑大孔中的晶体为球形，小孔中晶体的冰水平衡面为球面，公式（16）可写为：

Δ p = 2 γ_{s l} (\frac{1}{r_{2}} - \frac{1}{r_{1}})

（17）

式（17）与Everett给出的冰晶压力计算公式一致。从上述推导过程可以看出，Everett的公式是公式（10）的一个特例，只能在特定条件下使用，其使用条件为：1）大孔中存在已与孔壁接触的冰晶，且大孔中冰晶为球形；2）系统温度等于小孔中的冻结温度，小孔中冰晶与溶液处于平衡状态，小孔中溶液相对于大孔中的冰晶来说处于过冷状态；3）小孔中存在溶液供给大孔冰晶的生长。

2.4　胶囊状冰结晶压力

考虑一个球形晶体在圆柱形孔中生长，如图2所示。当球形晶体生长至与孔壁接触时，接触面上由于有结晶压力作用，该受荷面上的晶体生长受到结晶压力的抑制，而非受荷面上的晶体继续生长，晶体逐渐生长为端部为半球，侧边为圆柱的胶囊状晶体。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 胶囊状冰晶的生长示意图

Fig.2 Schematic diagram of crystal growth in a capsule pore

在圆柱形孔中的胶囊状晶体，端部球状部分（tip）与孔隙溶液平衡时（系统温度等于孔隙的冻结温度），其侧边圆柱部分（side）与孔壁接触，在该接触面上产生的结晶压力可由公式（16）计算：

Δ p = γ_{s l} (\frac{d A_{t i p}}{d V_{t i p}} - \frac{d A_{s i d e}}{d V_{s i d e}})

（18）

由于晶体端部为球形， $d A / d V = 2 / r$ ，侧部为圆柱有 $d A / d V = 1 / r$ ，代入式（18）可得：

Δ p = γ_{s l} \frac{1}{r}

（19）

公式（19）只适用于当系统温度等于孔隙的冻结温度时，圆柱形孔中胶囊状晶体与孔壁接触的侧边圆柱面产生的结晶压力计算。

上述分析表明，冰晶压力受溶液中水的活度、温度和冰晶表面曲率的共同影响，公式（10）是以晶体表面曲率为参量的表达式，公式（13）是以冻结温度为参量的表达式，二者在形式上有所区别但实质上是等价的。公式（8）、（17）及（19）均为公式（10）的特例，只有系统温度和晶体形状满足一定要求时才能使用。

3 晶体生长模型

3.1　晶体生长与孔壁的相互作用

结晶力计算公式得到的结晶压力为晶体生长过程中在晶体与孔壁接触面能够产生的最大结晶压力，该压力作用下受荷面上的晶体生长停止。指出对于弱固结的土体，由于本身强度较低，其内部无法达到和维持最大结晶压力，在孔隙间的结晶压力未达到最大结晶压力时，孔壁被晶体推动，土体即发生变形，晶体在受荷面继续生长。所以土体孔隙中晶体的结晶压力小于理论上的最大结晶压力，直接将理论上的最大结晶压力代入有效应力公式计算多孔材料变形会导致计算结果偏大。本模型中利用结晶压力来计算结晶力，利用结晶力和孔隙压力的关系来判断孔隙中晶体能否产生纵向生长。当多孔介质中的孔隙溶液达到某一临界状态时，开始有冰晶析出，析出的晶体不断生长，当晶体与孔壁接触后，在孔壁接触位置产生结晶压力且不断增大，当结晶力大于孔隙压力时，导致孔隙产生变形。基于该过程，本文将晶体生长过程划分为如下几个阶段：

（1）成核阶段：温度降低，当溶液温度达到临界过冷温度时，溶液中开始有冰晶成核；

（2）球形生长阶段［图3（a）］：温度持续降低，过冷度继续增大，冰晶体不断生长，为方便模型的计算，这里简化晶体为球形；

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 晶体生长示意图

Fig.3 Schematic diagram of crystal growth in a pore

（3）横向生长阶段［图3（b）］：当晶体与孔壁接触后，孔壁开始有结晶力作用，由于孔壁压力的约束作用，晶体只有横向生长，无纵向生长；

（4）横、纵向生长阶段［图3（c）］：当结晶力大于孔壁压力时，晶体开始出现纵向生长，此时晶体同时进行横向与纵向生长。当溶液过冷度小于0或者水完全转变为冰时，冰晶停止生长。

（5）接触生长阶段（图4）：当孔内晶核密度较大时，相邻的两个球形晶体会先彼此接触再与孔壁接触。让晶体与孔壁接触后并继续生长，也会出现两个胶囊形状的晶体彼此接触的情况。当晶体相互接触后，会融合为一个圆柱形晶体继续生长。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 晶体接触生长示意图

Fig.4 Diagram of crystal contact growth

3.2　过冷度与冰晶生长速率

降温过程中冰晶从溶液中析出的驱动力为过冷度，由过冷度计算冰晶成核率I，参考文献［9］中给出了具体的计算过程，这里不再赘述。假设冰晶晶核为球形，可由冰晶成核率计算冰晶生长速率mass_ice（kg·s^-1）：

\frac{1}{V} \frac{d m a s s_{i c e}}{d t} = \frac{4}{3} π {r^{*}}^{3} I ρ_{i c e}

（20）

式中：r*为晶核的临界半径；ρ_ice为冰晶密度。

3.3　孔隙中的晶体生长

在晶体与孔壁接触之前，晶体的生长为球形生长阶段，该过程中新的晶体半径可由新增结晶质量△m_s计算得到：

r_{s}^{t_{i + 1}} = {}^{3}{\sqrt[]{\frac{3 Δ m_{s}^{t_{i + 1}}}{4 π ρ_{c}} + (r_{s}^{t_{i}})^{3}}}

（21）

式中： $Δ m_{s}^{t_{i + 1}} = m_{s}^{t_{i + 1}} - m_{s}^{t_{i}}$ ； $ρ_{c}$ 为晶体密度。

当晶体与孔壁接触时，进入横向生长阶段，该阶段开始有结晶力产生F_s。本文采用公式（13）计算结晶压力 $Δ$ P，晶体与孔壁的接触面积A_s与孔壁所受结晶力F_s的计算公式为

F_{s} = Δ P \times A_{s}

（22）

A_{s}^{t_{i}} = 2 π r_{p}^{t_{i}} d_{s}^{t_{i}}

（23）

d_{s}^{t_{i}} = \frac{Δ m_{s}^{t_{i}} / ρ_{c} - 4 / 3 π (r_{p}^{t_{i}})^{3}}{π (r_{p}^{t_{i}})^{2}}

（24）

其中 $d_{s}^{t_{i}}$ 为晶体与孔壁接触长度，如图4所示， $r_{p}^{t_{i}}$ 为孔隙半径。

在横、纵向生长阶段孔隙半径r_p发生变化，开始出现孔隙变形，该阶段孔隙半径r_p等于晶体半径r_s，可由下式进行计算：

Δ m_{s}^{t_{i + 1}} / ρ_{c} = π \cdot d_{s}^{t_{i + 1}} [(r_{s}^{t_{i + 1}})^{2} - (r_{s}^{t_{i}})^{2}] + 4 / 3 π [(r_{s}^{t_{i + 1}})^{3} - (r_{s}^{t_{i}})^{3}]

（25）

3.4　数值模拟程序

3.4.1　模拟程序设定

本文对单个孔隙中的晶体生长过程进行模拟计算，对模型做如下设定：

（1）本模型模拟冻土的降温过程，孔隙中冰晶变形模拟程序温度变化如图5所示。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 模拟程序温度变化

Fig.5 Temperature of simulate program

（2）认为土体内部水分、晶核等均匀分布。晶核为球形，晶核半径等于临界半径为1 nm，晶核产生的数量正比于孔隙表面积，且平均分布于孔隙中，孔隙表面单位面积上的晶核数量为晶核密度。晶核的临界半径是晶核的生成时的最小尺寸，模型计算时假定为晶核的初始大小，且小于该尺寸的孔径中无晶核产生。

（3）孔隙及晶核密度设定：孔隙设定为圆柱形封闭孔，截取孔隙中只有一个晶核时的孔隙长度为一个孔单元，不同孔径孔隙的孔单元表面积相等。晶核密度设置的孔隙长度控制，孔隙长度越大对应的晶核密度越小，反之，孔隙长度越小对应的晶核密度越大。

（5）晶体横向时，晶体半径保持不变，晶体与孔壁接触面积增大。

（6）晶体的纵向生长只有在结晶力大于孔壁荷载的情况下产生，晶体纵向生长过程中，晶体与孔壁接触面积保持不变，孔隙半径增大。

（7）晶体的横、纵向生长阶段设定为晶体横向生长与纵向生长交替进行，进行纵向生长前需进行条件判断。

（8）土粒相对密度为2.71，土体孔隙比为0.34，土体干密度为1.75 g·cm^-3。

3.4.2　模拟计算流程

计算模拟程序采用MATLAB编写，程序流程如图6所示。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 冰晶生长模拟程序流程图

Fig.6 Flow chart of ice crystal growth simulation program

4 降温过程中孔隙变形分析

4.1　冰晶生长与孔隙变形分析

多孔介质中分布有大量不同孔径的孔，当其发生结晶变形时，每个孔内产生的晶体体积及每个孔的变形都不同，但当孔溶液的体积或浓度发生变化而导致孔内结晶体积发生变化时，虽然每个孔内结晶体积变化的大小不同，但其变化的规律一定是相同的，如随着孔溶液浓度的增加，每个孔内的结晶体积一定都是增加的，虽然每个孔内结晶体积的变化值不同。类似地，对于多孔介质中已发生冻胀变形的孔，当每个孔中结晶体积都增大时，每个孔的冻胀变形一定是增大的，虽然每个孔的变形大小不同，而此时，多孔介质的整体冻胀变形也必然是增大的，这说明单个孔隙的结晶变形与多孔介质的整体结晶变形虽然在数值上存在差异，但当环境参数发生变化时，二者的变化规律是相同的，这里对单个孔隙中的冰结晶变形过程进行分析。

考虑晶核密度为定值时，孔隙冻胀变形过程决定于孔隙中冰晶的短轴方向的生长过程，当冰晶的短轴长度大于孔隙直径时，孔隙开始发生变形，直到冰晶短轴方向生长结束，孔隙变形随之停止。孔隙变形量s可由孔隙中冰晶短轴长度l_s，和起始孔隙直径d₀计算得到（图7）：

s = l_{s} - d_{0}

（26）

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 结晶划分示意图

Fig.7 Diagram of crystal classification

当不考虑荷载作用时，根据几何关系，冰晶体积V_s与短轴长度l_s和长、短轴长度比（长宽比）λ之间有如下关系：

V_{s} = π {(\frac{l_{s}}{2})}^{3} (λ - \frac{2}{3})

（27）

将公式（26）代入公式（27）可得：

s = 2 \sqrt[3]{\frac{V_{s}}{π (λ - \frac{2}{3})}} - d_{0}

（28）

式（28）可以看到，孔隙的冻胀变形s同时受冰晶体积V_s、起始孔隙直径d₀和长宽比λ的影响，其中与冰晶体积V_s呈正相关关系，与起始孔隙直径d₀和长宽比λ呈负相关关系。从宏观上说，冻胀变形的大小与冰晶体积、孔隙大小和冰晶形状有关。

根据冰晶体积与起始孔隙直径的关系，将冰晶体积分成两部分，第一部分为填充结晶体积V_t，即只填充于孔隙中的晶体体积，不产生孔隙变形；第二部分为膨胀结晶V_p，这部分结晶为短轴长度超过原孔隙直径部分的晶体，该部分晶体的产生直接导致孔隙变形。定义在孔隙开始发生变形时刻，孔隙的冰晶体积为起胀结晶V_q，只有当冰晶体积大于起胀冰晶体积时，冻胀才会发生。从填充结晶与膨胀结晶关系的角度看，孔隙结晶中填充结晶越少，膨胀结晶越多，则孔隙的结晶变形越大，即晶体在孔隙中主要起填充效应还是膨胀效应，决定了孔隙的结晶变形的大小。这里定义孔隙中晶体体积与孔隙体积的比值为晶体填充率，结晶后孔隙增加的体积与原孔隙体积的比值为孔隙变形率，这两个参数能够分别反映晶体在孔隙中填充和膨胀效应的大小。显然，当结晶体积V_s一定时，膨胀结晶V_p所占比例越大，孔隙变形越大，反之，孔隙变形越小。起始孔隙直径d₀和长宽比λ对孔隙变形影响的机理，由图8可见，随着起始孔隙直径d₀和长宽比λ的增大，膨胀结晶V_p所占比例不断减小，从而导致孔隙变形减小。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 膨胀结晶比例与起始孔隙直径d₀和长宽比λ的关系曲线

Fig.8 Relationship between heaving crystals ratio and initial pore diameter d₀ and length-width ratio λ

由几何关系可得到起胀结晶V_q的计算公式为：

V_{q} = π {(\frac{d_{0}}{2})}^{3} (λ - \frac{2}{3})

（29）

由式（29）可以看到，起胀结晶体积V_q与起始孔隙直径d₀和冰晶长宽比λ呈正相关关系，在起始孔隙直径d₀一定，冰晶长宽比λ一定的情况下，起胀结晶体积V_q为定值。

由上述分析可知，冻胀变形的大小受冰晶体积、孔隙大小和冰晶形状的影响，当冰晶体积大于起胀结晶V_q时，孔隙开始膨胀，之后孔隙变形随着冰晶体积的增大而增大。而孔隙大小和冰晶形状影响了冰晶中膨胀结晶所占的比例，当孔隙直径越大，晶体长宽比越大时，膨胀结晶V_p所占比例越小，孔隙变形减小。

4.2　晶核密度对孔隙冰晶变形的影响

晶核密度是孔隙表面单位面积上的晶核数量，对材料内部的结晶应力及结晶变形有重要影响，由于晶核密度受环境温度，孔隙大小及分布、衬底的表面性质和晶体种类等多种因素的影响，缺少直接的测试手段，鲜有文献对晶核密度对孔隙变形的影机理响这一问题进行探讨。这里利用晶体生长模型，对晶体生长过程及晶核密度对孔隙变形的影响机理进行分析。

当晶核密度较大，单位面积孔隙上有较多的晶核同时生长，此时晶体的生长过程更趋向于先填充孔隙，再与孔壁发生作用导致孔隙变形。当晶核密度较小，单位面积孔隙上有较少的晶核同时生长，当晶体在孔壁处产生结晶压力时，孔隙中有大量空间尚未被填充。晶核密度影响了晶体在孔隙中的生长模式。

假设晶体的生长模式为先完全填充孔隙，然后在孔壁处产生结晶压力导致孔隙变形，该模式为完全填充模式。若晶体在孔壁处产生结晶压力时，孔隙中只有部分空间被填充，该模式为部分填充模式，如图9所示。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 孔隙中晶体生长模式示意图

Fig.9 Schematic diagram of crystal growth patterns in pores

许多学者以完全填充模式为基础，给出了结晶变形的计算方法^［14-16］。而实际上，完全填充模式只是晶体生长的一种特殊模式，该模式下晶体在孔隙中的填充率为100%，此时水相变成冰，体积膨胀8.8%。假设晶体生长模式为完全填充模式，当孔隙饱和度小于等于91.2%时，因为孔隙空间可以完全容纳冰晶的体积，理论上此时不发生冻胀。有大量试验数据表明，孔隙饱和度小于等于91.2%时，冻胀依然会发生，当孔隙体积明显大于析出的结晶体积时，会发生显著冻胀变形^［17-18］。完全填充模式也不能很好地解释多次冻融过程中冻胀变形的累加性。这些试验现象说明多孔材料中的晶体生长不能简单地用完全填充模式来描述，在某些条件下多孔材料中的晶体生长为部分填充模式，这时晶核密度的影响是需要考虑的。

选取孔径为2 μm，孔隙饱和度为100%，荷载为0.5 MPa，改变模型的晶核密度，得到不同晶核密度情况下，孔隙的结晶变形过程，如图10所示。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 孔隙结晶变形与晶核密度关系曲线

Fig.10 Relationship between pore crystal deformation and nucleus density

可以看到当晶核密度大于特征晶核密度X₁^*时，孔隙变形稳定为8.8%，从图10可知此时填充率为100%，此时冰晶在孔隙中为完全填充。此时的孔隙变形就是许多文献中所提到的，当材料内部水完全冻结时，冰晶的体积相对于水的体积膨胀8.8%。纯水溶液的冻结变形试验的结果与这一计算值是相符的。

当晶核密度小于特征晶核密度X₁^*，晶核密度越小，孔隙变形越大，此时晶体的生长模式为部分填充模式。由图10可知，此时随着晶核密度的减小，填充率不断减小，该阶段晶体的填充效应减弱，这是孔隙变形不断增大的原因。当晶核密度继续越小，孔隙变形不会无限增大，受荷载的抑制作用，当晶核密度小于特征晶核密度X₂^*时，孔隙结晶的填充率趋于稳定，所以此时孔隙变形也保持稳定。

4.3　荷载对孔隙冰晶变形的影响

有学者将荷载对冻胀变形的抑制作用，等效为荷载对冻土的压缩作用^［18］。当把问题细化到单个孔隙中，就会发现荷载对孔隙结晶变形的抑制作用无法等效为荷载对冰晶的压缩作用。由于在0 ℃时，冰的体积模量可以达到GPa级别^［19］，理论计算中冰被看作是无法压缩的。所以要弄清荷载对冻胀变形的抑制作用的机理，首先需要从单个孔隙入手进行分析。

选取孔径为2 μm，孔隙饱和度为100%，晶核密度为8 000 mm^-2，得到不同荷载作用下的孔隙的结晶变形，如图11所示。

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 孔隙结晶变形与荷载关系曲线图

Fig.11 Relationship between crystallization deformation and load in a pore

可以看到当荷载大于特征荷载P₁^*时，孔隙变形稳定为8.8%，此时填充率为100%，此时冰晶在孔隙中为完全填充。由于冰的体积模量巨大，理论计算上无法被压缩。由公式（13）可知，如果要抑制晶体的生成，所需很大压强，在-16 ℃时受荷面上冰晶应力达到20 MPa实际工程中如此高压状态几乎不存在，所以可以认为荷载对完全填充时的孔隙膨胀无抑制作用。

当晶核密度小于特征荷载P₁^*时，荷载越小，孔隙变形越大，此时晶体的生长模式为部分填充模式，与此同时，随着荷载的减小，填充率不断减小，该阶段荷载的抑制作用减弱，孔隙变形不断增大。当荷载小于特征荷载P₂^*时，荷载对填充率的影响很小，填充率开始趋于稳定，此时结晶变形也保持稳定。

由上述分析可知，荷载影响了结晶在孔隙中的填充率，这是由于荷载影响了孔隙中晶体的形状。如图12所示，荷载增加会导致晶体的长宽比增大，此时虽然冰晶的体积相同，但冰晶产生更多的横向生长，从而抑制了孔隙的结晶变形。

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 晶体长宽比与荷载关系曲线

Fig.12 Relationship between length-width ratio and load

4.4　孔径对孔隙冰晶变形的影响

许多研究认为当多孔材料中的孔隙空间越大，可以容纳越多的析出晶体，所以当结晶体积一定时，孔隙率越大，结晶变形越小。也有试验通过改变砂土级配控制土体的孔隙率，继而进行冻胀试验，试验结果表明其孔隙率与结晶变形之间没有明显联系^［20］。上述试验结果表明，现阶段对孔隙率影响结晶变形的机理的认知是不够全面和准确的。而孔隙率作为是多孔材料孔隙结构的宏观指标之一，其与结晶变形的关系本质上反映的是孔隙结构与结晶变形的关系。多孔材料的孔隙结构是指孔隙的大小、分布和连通性，而孔隙率只能反映多孔材料总的孔隙体积的大小，并不能完全反映多孔材料的孔隙结构。而孔隙结构对结晶变形有直接影响，比如有学者指出，结晶破坏与孔径大小直接相关，在大孔隙中，结晶破坏不容易产生，而在孔径为0.1~10 μm的孔隙中，更容易产生结晶破坏^［21］。由于孔隙结构不能简单用孔隙比来描述，应该针对孔隙的大小、分布及连通性对结晶变形的影响进行系统的分析。这一节主要利用晶体生长模型，分析孔径大小对孔隙结晶变形的影响。

选取晶核密度为5 300 mm^-2，孔隙饱和度为100%，荷载为0.2 MPa，改变模型的孔径，得到不同孔径的孔隙中结晶变形，如图13所示。

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13 孔隙结晶变形与孔径关系曲线

Fig.13 Relationship between crystallization deformation and pore diameter

当孔径大于特征孔径 $r_{1}^{*}$ 时，孔隙变形稳定为8.8%，此时冰晶在孔隙中为完全填充。由于20 ℃的纯水溶液完全冻结时，体积膨胀8.8%，所以该冻结变形试验可以看作发生在一个无限大饱和孔隙中，这说明一个无限大饱和孔隙中的冰晶也为完全填充，这与模拟的结果是吻合的。

当孔径小于特征孔径 $r_{1}^{*}$ 时，随着孔径不断减小，晶体的填充效应减弱，膨胀效应增强，孔隙变形快速增大。当孔径小于特征孔径 $r_{2}^{*}$ 时，受荷载的抑制作用，孔隙变形增速大幅降低，随着孔径不断减小，孔隙变形趋于稳定。

从图14可以明显看出，随着荷载的增大，特征孔径 $r_{2}^{*}$ 不断增大，且特征孔径 $r_{2}^{*}$ 对应的孔隙变形减小。荷载可以抑制小于特征孔径 $r_{2}^{*}$ 的孔隙的变形，而对于大于特征孔径 $r_{2}^{*}$ 的孔隙的变形，没有抑制作用。

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14 不同荷载条件下孔隙结晶变形与孔径关系曲线

Fig.14 Relationship between crystallization deformation and pore diameter under different loading conditions

5 结论

本文基于晶体热力学，对孔隙中冰晶对孔壁产生的最大结晶压力进行了研究，通过理论分析分别给出了不同形状晶体的结晶压力计算模型，并分析了经典结晶压力计算公式的使用条件。然后基于结晶动力学理论，考虑了晶体与孔壁的相互作用，建立了降温过程中孔隙冰晶生长模型，实现冰晶生长过程中的孔隙变形计算，分析了晶核密度、孔径大小、荷载和冰晶体积对孔隙结晶变形的影响机制，得到如下结论：

（1）起始孔隙直径和长宽比的增大对结晶变形抑制作用的机理在于减少了冰晶体积中膨胀结晶的比例。

（2）荷载对孔隙变形的抑制机制在于，荷载的增大迫使冰晶更多地横向生长（长宽比增大），导致膨胀结晶所占比例减小。

（3）孔隙中的晶体生长有完全填充模式和部分填充模式，在完全填充模式下，晶体在孔隙中的填充率为定值100%，在部分填充模式下，晶核密度、荷载和孔径的增大都会导致晶体在孔隙中的填充率增大，从而对孔隙结晶变形产生影响。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

， Lai

， Zhang

Thermo-hydro-salt-mechanical coupled model for saturated porous media based on crystallization kinetics

［J］. Cold Regions Science & Technology， 2016， 133： 94-107.