Page 89 - 水利学报2021年第52卷第5期

变饱和多孔介质内的水分运动可用水头型 Richards 方程描述，添加冰相的方程如下 ［9，13，27］ ：
C ∂h = Ñ × (kÑ∇ (h + i ) ) - ρ ∂θ i （4）
i
∂t ρ ∂t
w
式中：C 为比水容量；h 为基质势，m；k 为土体渗透系数，m/s；i 为重力项。
采用 van Genuchten 模型来描述未冻水含量与基质势、渗透系数的关系，方程如下 ［9，27］ ：
)
C = αm × (θ - θ ×Se 1 m × ( 1 - Se 1 m ) m （5）
1 - m s r
é
)
k (Se = k × Se × ê 1 - ( 1 - Se 1 m m ú ) ù 2 （6）
ë û
s
θ - θ 1/(1 - m ) m
(
Se = w r = 1 + |αh | ) （7）
θ - θ
s r
式中： α、m 为试验拟合参数； θ 、θ 分别为饱和、残余含水量；Se 为等效饱和度；k（Se）、 k 分别
s r s
为非饱和土、饱和土的渗透系数，m/s。
引入冰阻抗系数 I 来近似估算因冰的存在而引起的冻土渗透系数 k 的降低，方程如下 ［25-26］ ：
) θ ) -10θ i
k = k (Se × I ( ) = k (Se × 10 （8）
i
因颗粒表面能作用，冻土中的未冻水含量始终与温度保持动态平衡，方程如下 ［28］ ：
W = a|T | b （9）
u
式中： W 为未冻水的质量含水量；a、b 为实验参数；T 为含水量为 W 时对应的冻结温度（T），℃。
u u f
2.1.2 应力-应变控制方程 基于增量弹塑性理论，冻土的应力-应变方程如下 ［18］ ：
æ
}
{Dσ = [ ] {Dε - Dε p - v } ö （10）
D
} { } { Dε
è ø
{ }
式中： [ ] 为弹性矩阵； { Dε v } 为冻胀应变增量向量； Dε p 为塑性应变增量向量，由下式确定 ［18］ ：
D
{ } = λ p ∂{ } （11）
∂Q
p
Dε
σ
式中： λ 为塑性乘子；Q 为塑性势函数，满足相关联流动法则，等于屈服面函数 F，取 M-C 准则。
p
［9］
考虑冻土垂直和平行温度梯度的横观各向同性冻胀特征 ，并依据温度梯度方向实时修正主冻胀
［2］
方向 ，坐标转换至整体坐标系，冻胀应变向量方程如下：
2
ìε v ü é ê ê m 2 1 - ξ + n ξ ù ú ú
ï xx ï 2
ï
ï
ε υ ε v ý = (θ + θ - n ) ê ê 2 1 - ξ ú ú （12）
{ } = í yy w i 0 ê ê n + m 2
ï v ï 2 ú ú ξ
ï
ï
ε
î xy þ ê ê mn(1 - 3ξ ) ú ú
ë û
式中： ξ 为冻胀方向的分配权重，取 0.9；m、n 为温度梯度向量的方向余弦； n 为冻土初始孔隙率。
0
联 立 式（1）（4）（9）（10）构 成 渠 基 冻 土 的 温 度 场 、 水 分 场 和 变 形 场 的 耦 合 控 制 方 程 。 由 式（1）
（4）（9）可计算寒区渠道在外界环境作用下的温度场和冰、水含量分布，结合式（12）得到冻胀应变
参 数 ， 采 用 式（10）便 可 求 解 得 到 渠 道 的 应 力 变 形 分 布 ， 从 而 实 现 渠 基 冻 土 的 水 -热 -力 耦 合 冻 胀
分析。
2.2 渠道表面辐射换热模型
2.2.1 太阳辐射模型 太阳辐射具有明显的时空效应，其位置由赤纬 δ（°）和时角 ω（°）组成的赤道
坐标系及由太阳高度角 α 和方位角 γ 组成的地平坐标系决定，各参数计算公式如下 ［6，29］ ：
s
s
∘ ) ] （13）
δ = 23.45sin[ 360 × (284 + n /365
式中：n 为日序数，其中 1 月 1 日为 1，12 月 31 日为 365，表征不同季节日期的参数。
ω = 15(t - 12 ) （14）

— 591 —