黏性土比表面积较大,土颗粒吸附的强结合水含量占比较大,因而冻结过程中水的相变会持续产生^[19-20].

Tsytovich^[21]将冻结过程中的水相变划分为3个区间:剧烈相变区、过渡区和冻透区. 因颗粒均匀度高,黏土的剧烈相变区和过渡区温度区间持续时间相对较长. 某正冻黏土未冻水含量随负温变化如图1. 其中, θ_f为开始冻结温度; θ_n为进入冻透区温度; w₀和w_n为其相应的未冻水质量分数.负温θ_f时,黏性土进入冻结阶段,即开始产生剧烈相变; 负温θ_n时,土体冻透,土体中水冰等各相成分基本恒定,其热物理性质也趋于稳定^[22].

图1 正冻黏土中未冻水含量随负温变化
Fig.1 Variation of unfrozen water content in freezing clay soil with negative temperature

计算无热源的正冻黏土瞬态温度场时,若不考虑土中水分迁移,任意由光滑闭曲面S所围成的区域Ω单位时间内的热量变化Q₁, 应等于单位时间内通过曲面S流入(或流出)Ω内的热量Q₂, 即

Q₁=Q₂(7)

在单位时间内区域Ω的热量变化Q₁包括引起正冻黏土温度变化所需的热量和土中水结冰相变释放的潜热,即

Q₁=_Ω[(C_nf+C_L)ρ(θ)/(t)]dV(8)

Q₁=_Sλ_nf(gradθ)dS=

_Ωdiv(λ_nfgradθ)dV(9)

C_L=(ρ_w)/ρL(10)

其中,C_nf和λ_nf分别表示在区间[θ_n ～θ_f]上正冻黏土的比热和导热系数; C_L定义为冻土的相变热容,表示单位质量土温度每降低(或升高)1 ℃,土中未冻水冻结(或冰融化)释放(或吸收)的热量; ρ为黏土中水的密度; ρ_w为正冻黏土中水的密度.

由式(8)、式(9)和式(10)得

(C_nf+L(ρ_w)/ρ)ρ(θ)/(t)=div(λ_nfgradθ)(11)

由于正冻黏土中未冻水质量分数与负温保持着动态平衡的关系,即

w_u=w_u(θ), θ_n ≤θ≤θ_f(12)

通过未冻水含量测试试验拟合得到正冻黏土的未冻水含量关于负温θ在区间[θ_n～θ_f]上的函数w_u(θ), 可以推导出正冻黏性土的相变潜热在区间[θ_n～θ_f]上随未冻水含量变化的函数为

L(θ)=(335)/(1+w₀)w_u'(θ)Δθ, θ_n ≤θ≤θ_f(13)

其中, L(θ)为单位质量土体在负温θ时的相变潜热; w_u'(θ)为未冻水质量分数w_u(θ)的导数; Δθ为温度的微分.

式(13)中,左侧为状态变量,右侧是过程变量,故将区间[θ_n～θ_f]分为m等份, m值越大,则在θ_n≤θ_x≤θ_f上, L(θ_x)越接近L(θ), 于是有

L(θ_x)=

{(-335(θ_n-θ_f))/(2m+2mw₀)w_u'(θ_x), x=0,m

(-335(θ_n-θ_f))/(m+w₀)w_u'(θ_x), x=1,2,…,m-1(14)

正冻黏土中产生的相变潜热变成了若干个温度的积分,则在θ_n≤θ_x≤θ_f时,相应负温下的有效比热容C和有效导热系数λ为

C=(C_u+C_f)/2+(nρ_wL(θ_x))/(ρ(θ_n-θ_f))(15)

λ=λ_u,f=(λ_u+λ_f)/2(16)

其中,下角标u和f分别表示未冻区及冻实区, C_u,f和λ_u,f分别代表未冻区和冻实区的比热和导热系数.

将式(12)～(16)代入温度场计算方程(11),即可得到正冻黏土的瞬态温度场计算方程.

根据人工冻结试验的基本要求,设计模型尺寸为2.0 m×1.0 m×1.2 m. 依据土体保温要求,在模型槽内外壁均附设50 mm厚的保温板. 两根冻结管直径均为42 mm,冻结管水平贯穿模型槽短边方向,冻结管形心距模型箱外壁均为0.9 m,距离模型箱顶面及底面均为0.6 m,冻结管水平间距0.2 m,如图3.

图3 模型槽冻结管布置图
Fig.3 Freezing pipe layout of model slot

测温点均布置在以两根冻结管为法线的平面上,为减小模型箱体散热的影响,测点所在平面过冻结管中点,如图4. 其中,数字1～11代表测温点, L表示距离是5 cm, 2L和3L则分别表示距离是10 cm和15 cm.

图4 测温点所在平面布置图

模型试验用土为取自天津地铁2号线机场延长线上的粉质黏土,土密度为1.86×10³ kg/m³, w(H₂O)=26.8%. 实测正冻黏土中未冻水含量如表1. 拟合未冻水质量分数与温度的曲线如图5,土体材料和保温材料热物理参数见表2和表3,表4为基于实测数据计算所得的负温土体的潜热和有效比热容.

表1 实测正冻黏土的未冻水含量
Table 1 Measured content of unfrozen water in freezing clay

图5 未冻水质量分数随负温变化拟合曲线
Fig.5 Fitting curve of unfrozen water content with negative temperature

表2 土体材料热参数
Table 2 Thermal parameters of soil material

表3 保温材料热参数
Table 3 Thermal parameters of thermal insulation materials

假设冻结管与土体接触面为第1类边界条件,即土体与冻结管接触面的温度等于冷源温度. 设定冷源温度恒为-25 ℃,冻结管壁的实测温度如图6所示. 模型槽上表面为自然对流换热边界,对流换热系数为10 W/(m²·℃). 考虑到其余边界施加的保温材料并不能完全绝热,数值模拟过程中在模型槽其他面按照实际情况设置接触保温材料^[23]. 结合实测情况,设定保温材料外边界的温度为6.7 ℃,土体的初始温度为6.5 ℃.

表4 负温阶段黏性土中潜热计算
Table 4 Calculation of latent heat in cohesive soil at negative temperature stage

图6 冻结管外壁进出口温度随时间变化图
Fig.6 Temperature change of inlet and outlet of outer wall on freezing tube with time

图7所示的基于正冻黏土瞬态温度场数值计算方法(方法4),在土中水剧烈相变过后,方法4数值计算得到的温度曲线和实测曲线趋于完全重合. 在图7(b)及图7(i)中,在测点剧烈相变期间及之前阶段,模拟4得到的温度曲线和实测曲线仍存在少量误差,且表现为距冻结管越远,测点误差越大. 究其原因主要有:① 冷媒设定温度与实际供给温度存在差异,且土体各部分初始温度并不一致. ② 黏土未冻水含量测量存在误差. ③ 土中未冻水含量测试误差导致土体导热系数和比热产生持续误差,最终影响测试结果. 因此,距离冷源越远,温度变化越慢,土中未冻水含量测试精度对冻土的瞬态温度场计算影响越明显.

对比图7中4种模拟方法得到的各测点瞬态温度曲线可知,采用方法2和方法3计算得到的黏土瞬态冻结温度场远比方法1计算结果精确. 证实了相变潜热是影响土冻结过程的瞬态温度场计算准确性的重要因素. 方法4计算得到的黏土瞬态冻结温度场,比采用方法2和方法3更接近实测值. 因此,考虑黏土冻结过程中的未冻水含量变化,能够获取更加合理的土体瞬态冻结温度场.

方法1因未考虑土中水相变,其计算结果与实测值差异较大,且表现为距冻结管越远,模拟值与实测值差异越大,冻结管温度场叠加区域的实测值更接近模拟值. 采用方法2,即土中水在冰点全部相变,计算得到的测点瞬态温度曲线最不平滑,出现了个别突变点,计算结果收敛性降低,易出现相变遗漏. 相较之下,采用方法3,即等效比热容法计算得到的黏土瞬态冻结温度场,温度曲线更加平滑,且与实测值差异较小.

在方法3基础上,考虑了潜热在冻结过程中随土中未冻水含量的变化,使得方法4较其他计算方法获取的黏土瞬态冻结温度场更接近实测值. 因此,提出的黏土瞬态温度场计算方法不仅理论清晰,且计算精度更高.