人工冻土水热侵蚀问题室内试验研究

吴元昊; 孙旻; 洪泽群

doi:10.11779/CJGE2022S2038

人工冻土水热侵蚀问题室内试验研究 English Version

1.
中国建筑第八工程局有限公司, 上海 200135
2.
中国矿业大学力学与土木工程学院, 江苏徐州 221116

基金项目:

上海市青年科技英才扬帆计划 21YF1452300

详细信息

作者简介:
吴元昊(1991—)，男，河南驻马店人，博士，主要从事冻结法及隧道地下工程的研究工作。E-mail: wu_yuanhao163@163.com

中图分类号: TU43
计量
- 文章访问数: 70
- HTML全文浏览量: 13
- PDF下载量: 20
出版历程
- 收稿日期: 2022-11-30
- 网络出版日期: 2023-03-26
- 刊出日期: 2022-11-30

Laboratory tests on fluvial-thermal erosion of artificial frozen soil

1.
China Construction Eighth Engineering Division Co., Ltd., Shanghai 200135, China
2.
School of Mechanics and Civil Engineering, China University of Mining and Technology, Xuzhou 221116, China

摘要

摘要: 对人工冻结工程灾害过程中发生的人工冻土水热侵蚀问题进行了概念上的定义，并指出该问题的实质是双移动边界位置的确定。采用自行研制的冻土侵蚀室内模拟试验装置对该问题进行研究，介绍了试验装置及测点布置情况，给出双移动位置的定量判定标准。选择影响水热侵蚀的水流速率、冻土初始温度、土体含水率3个主要因素作为试验变量，采用控制变量法进行试验分组，得到不同情况下双边界随时间的移动变化情况。试验研究发现：双边界移动速率随流速呈现两阶段变化规律，即先快速增加后平稳趋势，表明冻土的存在使侵蚀速率大大减缓；冻土初始温度对于边界移动的影响并不明显，3种情况下侵蚀边界移动情况基本一致；土体饱和与否对水热侵蚀边界移动有很大影响，土体未饱和时对应的潜蚀速率明显大于饱和时的潜蚀速率。
- 冻结工程 /
- 灾变机理 /
- 水热侵蚀 /
- 双移动边界 /
- 试验研究
Abstract: The problem of fluvial-thermal erosion of artificial frozen soil during the process of ground freezing disasters is conceptually defined, and its essence is the determination of the position of the double-moving boundary. The self-developed indoor simulation test device for the erosion of frozen soil is used to study the problem. The test device and the arrangement of measuring points are introduced, and the quantitative judgment standard for the double–moving boundary is given. Three main factors affecting the fluvial erosion are selected: water flow velocity, initial temperature of frozen soil and moisture content of soil. The control variable method is used to group the tests. The movements of the double boundaries over time under different conditions are obtained. It is found that the movement rate of the double boundaries shows a two-stage change pattern with the flow velocity. The movement rate increases rapidly and then slows down, indicating that the existence of the frozen soil greatly slows down the erosion rate. Under three different temperature conditions, the movement of erosion boundary is basically the same, indicating that the initial temperature of the frozen soil has few influences on the boundary movement. Whether the soil is saturated or not has a great influence on the movement of erosion boundary, and the corresponding erosion rate when the soil is not saturated is obviously greater than that when it is saturated.
- ground freezing work /
- catastrophic mechanism /
- fluvial-thermal erosion /
- double-moving boundary /
- laboratory test

HTML全文

0. 引言

20世纪末以来,国内外陆续建设成一批超高面板坝工程,在取得成功建设经验的同时,也发现有不少工程在运行期发生了混凝土面板沿纵接缝的挤压破损问题^[1],例如,中国的天生桥一级面板坝（坝高178 m）^[2-3]、水布垭面板坝（坝高233 m）^[4]等,国外的默霍尔面板坝（Mohale,坝高145 m）、巴拉格兰德面板坝（Barra Grande,坝高185 m）和肯柏诺沃面板坝（Campos Novos,坝高202 m）等^[5-7]。

上述实例是超高面板堆石坝发生最为普遍的一种面板结构破损。文献[8]总结了这些工程中面板发生挤压破损的现象,发现这种挤压破损现象具有都发生在河谷中央面板的顶部部位、压性纵缝的两侧以及通常首先发生在面板的表层等特征。图1（a）,（b）分别给出了天生桥一级面板坝面板发生挤压破损的照片和位置图。

图 1 天生桥一级面板坝面板挤压破损及挤压应力计算值

Figure 1. Extrusion damages and calculated extrusion stresses of concrete face in Tianshengqiao-1 CFRD

下载: 全尺寸图片幻灯片

混凝土面板通常被看作是堆石体上的柔性薄板结构,其应力的大小和分布取决于堆石坝体变形的大小分布。在坝轴线方向,坝体堆石体发生由两岸指向河谷方向的位移,从而使河谷中部的面板发生轴向挤压,这被认为是面板发生挤压破损的主要原因。

但是,根据图1（c）所示的计算结果可以发现,计算结果和实际面板挤压破损之间存在诸多的矛盾：①实际挤压破损并没有发生在计算挤压应力最大的位置。有限元计算的最大坝轴向挤压应力一般总是位于河谷中央面板的中部或中下部,而面板的实际挤压破损却总是发生在河谷中央面板的顶部。②计算所得面板的最大挤压应力一般不超过15 MPa,尚远低于面板混凝土的极限抗压强度。即使考虑面板顶部相对不利的受力条件,面板轴向抗压也应是安全的。③天生桥一级面板堆石坝已安全运行20 a,坝体后期流变变形的量值已经不大,应该不再是造成面板轴向挤压的主要因素。

通过上面的讨论可以发现,对于面板的挤压破损现象,除了由于坝体变形所导致的面板挤压之外,还应该有其它重要的影响因素。

天生桥一级面板堆石坝最大坝高178 m,是中国首座超高面板堆石坝。该坝位于贵州和广西界河南盘江中下游,在北纬25度附近,夏季太阳辐射强度大。在面板堆石坝中,面板位于大坝上游面,其应力受外界环境温度变化的影响较大。在发生大幅升温时也可产生很大的坝轴向挤压应力。据统计,自2003年—2012年,天生桥一级面板堆石坝共发生8次面板挤压破损,均发生在夏初或夏中,其中发生在5月下旬2次、6月3次、7月3次。根据现场的观测资料,在太阳热辐射作用下,夏季天生桥面板坝面板的表面温度最高可达50℃～60℃。

通过以上分析可以发现,夏季太阳热辐射在面板中产生高温,是否会在面板中产生较高的温度应力,从而诱发面板发生挤压破损,是一个值得深入研究和探讨的问题。关于面板温度应力,以往的研究重点多集中在施工期及面板裂缝等问题,例如浇筑时的气温、入仓温度、水化热和气温变化的影响^[9-12],施工期太阳辐射的影响^[13]等。王子健等^[14]考虑运行期季节性气温变化,从结构应力和温度应力的角度分析了公伯峡面板堆石坝运行期面板产生裂缝的原因。此外,还有考虑温度对面板混凝土模量及堆石料变形影响的相关研究^[15-16]。目前,有关面板坝运行期尤其是太阳辐射温度应力的研究成果还较为少见。

本文基于对偶mortar有限元多体接触计算方法,发展了一种瞬态温度场分析与应力变形分析的热-力多物理场耦合计算方法。对一理想高面板坝进行有限元计算,分析了夏季太阳热辐射作用对混凝土面板应力的影响。

1. 热-力耦合对偶mortar有限元法

在面板坝有限元计算分析中,面板应力的准确计算是个难题。在面板坝中存在面板与垫层、面板与面板以及面板与周边缝等接触关系,因而是一个典型的复杂多体接触系统。采用传统的接触面单元法,面板应力计算的精度较差。此外,面板厚度很薄,造成面板单元划分的几何形态一般很差^[17]。本文要进行热-力耦合计算,需分析温度及其应力沿面板厚度方向的分布,这要求在面板厚度方向划分更多层的单元,从而会使得面板单元划分的矛盾更加突出。

为了解决上述的关键难题,本文发展了基于对偶mortar元的计算接触力学方法。相比传统的接触面单元法,该法对位移不连续现象的处理具有本质的优越性。此外,对每个接触物体还可以单独划分有限元网格。也即,采用这种接触方法进行面板坝计算时,坝体和面板的单元可分别单独划分,不要求两者之间协调。这样很容易实现对面板网格的精细划分,从而提高面板应力的计算精度。

mortar元^[18]作为最流行的一种面-面接触方法,已成功嵌入ABAQUS和ANSYS等商业软件,可用以应对接触求解精度的挑战,但由于矩阵性态的影响,大规模计算的求解效率仍需进一步考察。对偶mortar元^[19-20]可保持与mortar元同精度,且同时提升求解效率,为面板应力的精细计算提供了一种新途径。

1.1 基于对偶mortar元的非稳定温度场计算方法

对运行期面板堆石坝坝体和面板的非稳定温度场,采用瞬态热传导方程^[21-22]：

$c ρ \dot{θ} - k \nabla^{2} θ = 0$ ,

(1)

式中,c,ρ,θ和k为比热、密度、温度和热传导系数。符号 $\nabla^{2}$ 为对空间坐标的二阶偏导算子。三类边界条件为（图2）：

图 2 非稳定温度场的三类边界条件

Figure 2. Three boundary conditions in thermal analysis

下载: 全尺寸图片幻灯片

： ， ： ， ： ，

$\begin{array}{l} Γ_{θ} ： θ = \hat{θ} ， \\ Γ_{q} ： q \cdot n_{q} = {\hat{q}}_{n} ， \\ Γ_{r} ： q \cdot n_{r} = h_{c} (θ - θ_{a}) - α_{s} q_{s} ， \end{array}}$

(2)

式中, $\hat{θ}$ 为Dirichlet边界 $Γ_{θ}$ 上给定的温度。 $q$ 为边界热流。 $n_{q}$ , ${\hat{q}}_{n}$ 分别为Neumann边界 $Γ_{q}$ 的单位法向量和给定的热流。 $n_{r}$ 为Robin边界Γ_r的单位法向量,h_c,θ_a分别为对流系数和空气温度, $α_{s}$ 为表面吸热系数,q_s为与太阳辐射相对应的热流。

本文采用基于对偶mortar元的计算接触力学方法,因而还存在不同物体间的接触界面。在本文的计算中,为使面板具有较高计算精度而采用较密网格,垫层与堆石则采用较稀网格。对于三者之间的非协调网格接触界面Γ_C,需满足如下连续性条件：

$θ_{s} = θ_{m}$ ,

(3)

式中, $θ_{s}$ , $θ_{m}$ 为从面和主面的温度。

基于约束变分原理,使用Lagrange乘子法施加上述连续性条件,可得控制方程的弱形式为

$\begin{array}{l} \int_{Ω^{t}} [δ θ (c ρ \dot{θ}) + \nabla (δ θ) \cdot (k \nabla θ)] d Ω + \\ \int_{Γ_{r}} δ θ (h_{c} θ) d Γ + \int_{Γ_{C}} (δ θ_{s} - δ θ_{m}) λ d Γ = \\ \int_{Γ_{r}} δ θ (h_{c} θ_{a} + α_{s} q_{s}) d Γ - \int_{Γ_{q}} δ θ {\hat{q}}_{n} d Γ \end{array}$ ,

(4)

式中,Ω^t为当前计算区域。

相应地,接触界面连续性条件的弱形式为

$\int_{Γ_{C}} δ λ (θ_{s} - θ_{m}) d Γ = 0$ 。

(5)

温度θ及其变分δθ采用有限元线性单元的形函数N进行插值,Lagrange乘子 $λ$ 及其变分 $δ λ$ 则采用如下的对偶基函数 $Ψ$ 进行插值^{[19-20, 23]}：

$\int_{Γ_{C}} Ψ_{j} N_{k} d Γ = δ_{j k} \int_{Γ_{C}} N_{k} d Γ$ ,

(6)

式中,δ_jk为Kronecker函数。

采用后向欧拉法进行时间插值,对式（4）进行时间积分,可得基于对偶mortar有限元法的非稳定温度场计算格式：

$\begin{matrix} H θ + L^{T} λ - C^{T} λ = Q ， \\ L θ_{S} - C θ_{M} = 0 ， \end{matrix}}$

(7)

式中,

$\begin{matrix} H = \int_{Ω^{t}} (c ρ N^{T} N + τ k \nabla N^{T} \nabla N) d Ω + \\ \int_{Γ_{r}} τ h_{c} N^{T} N d Γ \end{matrix}$ ,

(8a)

$\begin{array}{l} L = \int_{Γ_{C}} τ Ψ^{T} N_{S} d Γ ， C = \int_{Γ_{C}} τ Ψ^{T} N_{M} d Γ \\ Q = \int_{Ω^{t}} c ρ N^{T} N θ_{his} d Ω - \int_{Γ_{q}} τ N^{T} {\hat{q}}_{n} d Γ + \end{array}$ ,

(8b)

$\int_{Γ_{r}} (τ h_{c} N^{T} θ_{a} + τ α_{s} N^{T} q_{s}) d Γ$ ,

(8c)

其中,τ表示相邻时间步的时间间隔。

将结点分为3个子集：从点S、主点M以及其他点R。式（7）可等价地写为以下的子块形式：

$[\begin{matrix} H_{RR} & H_{RS} & H_{RM} & 0 \\ H_{SR} & H_{SS} & H_{SM} & L^{T} \\ H_{MR} & H_{MS} & H_{MM} & - C^{T} \\ 0 & L & - C & 0 \end{matrix}] {\begin{matrix} θ_{R} \\ θ_{S} \\ θ_{M} \\ λ \end{matrix}} = {\begin{matrix} Q_{R} \\ Q_{S} \\ Q_{M} \\ 0 \end{matrix}}$ 。

(9)

可以看到,左端项矩阵最后一行的对角元为零元素,导致矩阵失去正定性,一般需采用稳定双共轭梯度法（Bi-conjugate gradients stabilized method）等特殊方法进行求解,且求解效率也不高。

注意到式（6）给出的对偶基函数 $Ψ$ 将使得L矩阵成为严格对角矩阵,因此可使用矩阵变换方法^[23],将式（9）中的第2行和第4行消去,同时将第2列和第4列也消去,从而得到：

$[\begin{matrix} H_{RR} & {\tilde{H}}_{RM} \\ {\tilde{H}}_{RM}^{T} & {\tilde{H}}_{MM} \end{matrix}] {\begin{matrix} θ_{R} \\ θ_{M} \end{matrix}} = {\begin{matrix} Q_{R} \\ Q_{M} + G^{T} Q_{S} \end{matrix}}$ ,

(10)

其中,

$\begin{array}{l} \begin{array}{l} G = L^{- 1} C ， \\ {\tilde{H}}_{RM} = H_{RS} G + H_{RM} ， \end{array} \\ {\tilde{H}}_{MM} = H_{MS} G + H_{MM} + G^{T} H_{SM} + G^{T} H_{SS} G 。 \end{array}}$

(11)

在式（11）中,矩阵 $G$ 包含矩阵L的逆。如果将 $Ψ$ 直接取为温度的插值函数N（对应常规mortar元^[18]）,此时,矩阵L的求逆代价极高,往往不如直接求解式（9）。而在对偶mortar元中,由于L的严格对角特性,其求逆代价可忽略不计,因此可选择求解式（10）,此时不仅左端项矩阵恢复了正定性,而且矩阵规模也大为缩减,从而计算效率可显著提升。这一优势完全得益于式（6）中对偶基函数 $Ψ$ 的定义,这也正是对偶mortar元与常规mortar元的本质区别。

1.2 基于对偶mortar元的热-力耦合计算方法

对应力变形计算,控制方程为

$\nabla \cdot σ + b = 0$ ,

(12)

式中, $σ$ 为应力张量, $b$ 为体积力。

在计算中,存在面板-垫层以及面板-面板等之间的非线性摩擦接触。摩擦接触可描述为如下的法向和切向KKT条件^{[20, 23]}：

$\begin{array}{l} \begin{matrix} 法向： & p_{n} \leq 0 ， & g_{n} \geq 0 ， & p_{n} g_{n} = 0 \end{matrix} ， \\ \begin{matrix} 切向： & κ_{τ} \leq 0 ， & ι_{τ} \geq 0 ， & κ_{τ} ι_{τ} = 0 \end{matrix} ， \end{array}}$

(13)

式中, $p_{n}$ 为法向接触应力,以压为正。g_n为法向间隙,以脱开为正。κ_τ和ι_τ定义为

$\begin{array}{l} κ_{τ} = ‖ f_{τ} ‖ - μ | p_{n} | ， \\ {\bar{υ}}_{τ} + ι_{τ} f_{τ} = 0 ， \end{array}}$

(14)

式中, $f_{τ}$ 为切向摩擦应力,μ为摩擦系数, ${\bar{υ}}_{τ}$ 为切向相对运动速度。可见,KKT条件规定了法向应满足不可贯入条件和无拉应力条件,而切向应满足库伦定律中的黏结条件和摩擦运动学条件。

采用Lagrange乘子法施加KKT条件,可得控制方程的弱形式为

$\begin{matrix} \int_{Ω^{t}} 〈 \nabla (δ u), σ 〉 d Ω + \int_{Γ_{C}} (δ u_{s}^{T} - δ u_{m}^{T}) \tilde{R} \tilde{λ} d Γ \\ = - \int_{Ω^{t}} δ u \cdot b d Ω - \int_{Γ_{σ}} δ u \cdot \hat{t} d Γ \end{matrix}$ ,

(15)

其中, $u$ 为位移。 $u_{s}$ 和 $u_{m}$ 分别为从面和主面上的位移。 $\tilde{λ}$ 为应力变形计算中的Lagrange乘子。 $\hat{t}$ 为Neumann边界Γ_σ上的给定面力。矩阵 $\tilde{R}$ 表示由从面法向n和切向 $τ_{1}$ , $τ_{2}$ 构成的正交旋转矩阵：

$\tilde{R} = [τ_{1}, τ_{2}, - n]$ 。

(16)

相应地,KKT条件的弱形式为

$\begin{array}{l} \int_{Γ_{N}} δ {\tilde{λ}}^{T} {\tilde{R}}^{T} (u_{s} - u_{m}) d Γ = 0 ， \\ \int_{Γ_{H}} δ u_{s}^{T} {\tilde{δ}}_{τ} \tilde{λ} d Γ = \int_{Γ_{H}} δ u_{s}^{T} {\tilde{R}}_{τ}^{T} μ | {\tilde{λ}}_{n}^{tr} | τ^{tr} d Γ ， \\ \int_{Γ_{H}} δ {\tilde{λ}}^{T} {\tilde{R}}_{n}^{T} (u_{s} - u_{m}) d Γ = 0 ， \end{array}}$

(17)

式中,Γ_N,Γ_H分别为处于黏结状态和滑动状态的边界。 ${\tilde{λ}}_{n}^{tr}$ , $τ^{tr}$ 分别为预估的法向接触应力和从面切线方向。矩阵 ${\tilde{R}}_{τ}$ , ${\tilde{R}}_{n}$ 和 ${\tilde{δ}}_{τ}$ 为

$\begin{array}{l} {\tilde{R}}_{τ} = [τ_{1}, τ_{2}, 0] ， \\ {\tilde{R}}_{n} = [0, 0, - n] ， \\ {\tilde{δ}}_{τ} = {\tilde{R}}^{T} {\tilde{R}}_{τ} 。 \end{array}}$

(18)

使用Newton-Raphson方法处理材料非线性,则应力张量 $σ$ 可线性化为

$σ \approx \frac{\partial σ}{\partial ε} Δ ε + σ_{his}$ ,

(19)

式中, $ε$ 为应变张量, $σ_{his}$ 为历史应力。应变增量 $Δ ε$ 为位移增量 $Δ u$ 和温度增量Δθ的函数：

$Δ ε = \frac{\partial ε}{\partial u} Δ u - α_{θ} Δ θ I_{3 \times 3}$ ,

(20)

其中, $α_{θ}$ 为温度线膨胀系数, $I_{3 \times 3}$ 为二阶单位球张量。

位移u及其变分 $δ u$ 采用有限元线性单元的形函数 $\tilde{N}$ 进行插值,Lagrange乘子 $\tilde{λ}$ 及其变分 $δ \tilde{λ}$ 则采用类似式（6）的对偶基函数 $\tilde{Ψ}$ 进行插值：

${\tilde{Ψ}}_{j} = [Ψ_{j}, Ψ_{j}, Ψ_{j}]$ ,

(21)

式中, $Ψ_{j}$ 也就是式（6）所给出的对偶基函数, ${\tilde{Ψ}}_{j}$ 则表示此处用于插值结点乘子 ${\tilde{λ}}_{j}$ 的相应函数。

至此,可得基于对偶mortar有限元法的应力变形计算格式：

$\begin{array}{l} K Δ u + D^{T} R \tilde{λ} - M^{T} R \tilde{λ} = F ， \\ R_{NN}^{T} D_{NN} Δ u_{N} - R_{NN}^{T} M_{NM} Δ u_{M} = r_{N} ， \\ R_{n}^{T} D_{HH} Δ u_{H} - R_{n}^{T} M_{HM} Δ u_{M} = r_{H}^{n} ， \\ δ_{τ} D_{NN}^{T} {\tilde{λ}}_{H} = r_{H}^{τ} ， \end{array}}$

(22)

式中,下标N和H分别表示从面上处于黏结和滑动状态的相应结点集合,下标M表示主面上的结点集合。矩阵R,R_n和 $δ_{τ}$ 分别为矩阵 $\tilde{R}$ , ${\tilde{R}}_{n}$ 和 ${\tilde{δ}}_{τ}$ 的结点平均形式。矩阵K,D,M,F以及向量r分别为

$K = \int_{Ω^{t}} (\nabla {\tilde{N}}^{T} \frac{\partial σ}{\partial ε} \nabla \tilde{N}) d Ω$ ,

(23a)

$D = \int_{Γ_{C}} {\tilde{Ψ}}^{T} {\tilde{N}}_{S} d Γ ， M = \int_{Γ_{C}} {\tilde{Ψ}}^{T} {\tilde{N}}_{M} d Γ$ ,

(23b)

$\begin{matrix} F = \int_{Ω^{t}} (α_{θ} \nabla {\tilde{N}}^{T} \frac{\partial σ}{\partial ε} m N Δ θ) d Ω - \\ \int_{Ω^{t}} (\nabla {\tilde{N}}^{T} σ_{his} + {\tilde{N}}^{T} b) d Ω - \int_{Γ_{σ}} {\tilde{N}}^{T} \hat{t} d Γ \end{matrix}$ ,

(23c)

$r_{N} = - R_{NN}^{T} (D_{NN} x_{N} - M_{NM} x_{M})$ ,

(23d)

$r_{H}^{n} = - R_{n}^{T} (D_{HH} x_{H} - M_{HM} x_{M})$ ,

(23e)

$r_{H}^{τ} = R_{τ}^{T} D_{HH}^{T} e_{n}^{T} | {\tilde{λ}}_{H}^{tr} | \int_{Γ_{H}} μ τ^{tr} d Γ$ ,

(23f)

式中, $R_{τ}$ 为矩阵 ${\tilde{R}}_{τ}$ 的结点平均形式。m为二阶单位球张量 $I_{3 \times 3}$ 的向量表达,规则同应力张量的向量表达。向量x为结点坐标。

采用混合切向接触条件^[23]解决面板相互接触所引发的算法不稳定问题,采用矩阵变换方法^[23]解决鞍点问题,则式（22）可改写为

$\begin{array}{l} (P^{T} K P) y = P^{T} F - D^{T} \tilde{λ} ， \\ 滑动： e_{τ_{1}}^{T} {\tilde{λ}}_{j} = μ e_{n}^{T} | {\tilde{λ}}_{j}^{tr} | ， \\ 滑动： e_{τ_{2}}^{T} y_{j} = 0 ， \\ 滑动： e_{n}^{T} y_{j} = - n_{j}^{T} (x_{j} - D_{j j}^{- 1} \sum_{k} M_{j k} x_{k}) ， \\ 黏结： y_{j} = - R_{j j}^{T} (x_{j} - D_{j j}^{- 1} \sum_{k} M_{j k} x_{k}) 。 \end{array}}$

(24)

式中y为变换后的待求未知量, $e_{n}^{}$ ,e_τ为单位坐标基,P为变换矩阵,

$Δ u = P y$ ,

(25a)

$\begin{array}{l} e_{τ_{1}} = {[1, 0, 0]}^{T} ， \\ e_{τ_{2}} = {[0, 1, 0]}^{T} ， \\ e_{n} = {[0, 0, 1]}^{T} ， \end{array}}$

(25b)

$P = [\begin{matrix} I_{RR} & 0 & 0 & 0 \\ 0 & I_{MM} & 0 & 0 \\ 0 & D_{NN}^{- 1} M_{NM} & R_{NN} & 0 \\ 0 & D_{HH}^{- 1} M_{HM} & 0 & R_{HH} \end{matrix}]$ 。

(25c)

可见,经过变换之后,式（24）中的接触约束成为了完全解耦且正交于向量y或 $\tilde{λ}$ 的形式,因此可直接等同于施加Dirichlet或Neumann边界条件。

2. 热-力耦合算法流程及计算程序

（1）求解式（10）,得到温度场θ。

（2）将θ代入式（23c）,得到右端项F。

（3）求解式（24）,得到解向量y。

（4）根据式（25a）,计算得到位移场u。

（5）根据式（24）,反算得到接触应力 $\tilde{λ}$ 。

（6）根据式（19）,（20）,更新历史应力 $σ_{his}$ 。

（7）根据PDASS方法^[24],更新接触状态。

（8）若残差较小,则认为计算收敛,并进入下一时间步；否则,返回至第（1）步进行迭代计算。

根据上述计算方法,自主开发了相应的有限元数值模拟程序系统。

3. 200 m级理想面板堆石坝计算分析

3.1 计算模型与参数

参照古水面板坝的设计资料,考虑一坝高为200 m的理想坝,其坝顶长600 m,坝顶宽10.26 m,河谷底宽48 m,上游坝坡1∶1.5,下游坝坡分别为1∶1.7和1∶1.4。正常蓄水位和死水位分别取182.2,153.6 m。

材料分区如图3所示,包括主堆石料3B1和3B2,次堆石料3C,排水区堆石料3D,过渡区堆石料3A1和3A2,垫层堆石料2A。面板厚度沿高程线性变化,顶部取0.5 m,底部取1.1 m。面板纵缝间距布置为12 m,共将面板划分为50块。

图 3 理想坝的材料分区

Figure 3. Material zones and filling stages of ideal dam

下载: 全尺寸图片幻灯片

根据材料分区和施工分级,剖分了如图4所示的有限元计算网格,包括194万个节点,146万个单元,579万自由度。其中,为了能准确分析温度及应力沿面板厚度方向的分布,在面板厚度方向划分了6层的单元,相应地,对面板单元在平面方向也进行了加密处理,平面内的单元尺寸约为1 m左右。

图 4 三维计算网格

Figure 4. Three-dimensional computational mesh

下载: 全尺寸图片幻灯片

本文采用所发展的基于对偶mortar元的计算接触力学方法进行计算。为了降低计算规模,坝体部分单元划分尺寸仍取10 m左右的量级。另将垫层作为面板和堆石间的接触过渡区。这样,面板、垫层和坝体可根据各自的需要单独划分有限元网格,不要求计算网格的协调,很容易实现面板网格的精细划分。这也是计算接触力学方法的一大优势。图4右侧局部放大图给出了它们之间的接触过渡关系。这样,在本文的计算中,共包括有坝体、垫层和50块面板共52个物体,它们共同构成一个复杂的多体接触系统。

参照文献[14,22],选定了材料热力学性能相关的计算参数如表1所示。

表 1 材料热学性能参数

Table 1. Thermal parameters of different materials

材料	cρ/(10⁶J·m^-3·℃^-1)	k/(J·m·s^-1·℃^-1)	h_c/(J·m^-2·s^-1·℃^-1)	$α_{s}$	$α_{θ}$ /(10^-7℃)
主堆石	1.63	1.227	12	0.60	8.5
次堆石	1.63	1.227	12	0.60	8.5
排水区	1.63	1.227	12	0.60	8.5
过渡料	1.96	1.717	12	0.60	3.0
垫层料	1.99	1.472	12	0.60	3.0
面板	2.45	2.453	20	0.65	1.0

下载: 导出CSV

| 显示表格

本文进行了夏季典型一天的非稳定温度场计算。为简化计,假定6∶00时为初始条件,温度取20℃。对图2中所示的边界条件：对水面以下的面板,采用给定水温20℃；对大坝底部与坝基接触的部分,采用绝热边界；对暴露在阳光下的部分,即坝顶、下游堆石表面以及水面以上的面板,均采用热交换边界。

对热交换边界,考虑太阳热辐射和气温的变化。太阳辐射热流量 $q_{s}$ 及气温 $θ_{a}$ 均假设为余弦分布^[22]：

$\begin{array}{l} q_{s} = \frac{π Q_{s}}{2 t_{s}} \cos \frac{π {\bar{t}}_{mid}}{t_{s}} ， \\ θ_{a} = θ_{a0} + θ_{ap} \cos \frac{π {\bar{t}}_{mid}}{t_{s}} ， \end{array}}$

(26)

式中, $Q_{s}$ 为日均热量,取为3.28×10⁷ J/m²,对应北纬30°～40°的夏季值。t_s为日照持续时间,取为12 h。 ${\bar{t}}_{mid}$ 为当前时刻与日照最强时刻（12∶00）的差值。初始气温 $θ_{a0}$ 和气温变幅 $θ_{ap}$ 均取为20℃。

在应力变形计算中,堆石料采用邓肯张E-B模型及七参数流变模型^[25],具体参数由古水面板坝相关室内试验确定,分别见表2,3所示。混凝土采用线弹性模型,杨氏模量E = 3.0×10¹⁰ Pa,泊松比 $ν$ = 0.2。面板与垫层之间的摩擦系数取为0.75,面板与面板之间的摩擦系数取为0.7。

表 2 坝料的邓肯张E-B模型参数

Table 2. Parameters of E-B model for rockfill materials

坝料	φ/(°)	∆φ/(°)	K	n	R_f	K_b	m
主堆	55.5	11.3	1350	0.28	0.80	780	0.18
次堆	53.0	11.0	1000	0.26	0.79	700	0.16
排水	55.0	12.2	1300	0.31	0.79	800	0.12
过渡	53.5	10.7	1250	0.31	0.78	720	0.16
垫层	54.4	10.6	1200	0.30	0.75	680	0.15

下载: 导出CSV

| 显示表格

表 3 坝料的流变模型参数

Table 3. Material parameters of creep model for rockfills

坝料	α	b	c	d	m₁	m₂	m₃
主堆	0.0012	0.0008	0.96	0.0012	0.39	0.41	0.63
次堆	0.0017	0.0010	0.74	0.0015	0.47	0.48	0.66
排水	0.0017	0.0008	0.96	0.0012	0.39	0.41	0.63
过渡	0.0012	0.0008	0.97	0.0013	0.39	0.41	0.63
垫层	0.0012	0.0007	0.95	0.0011	0.40	0.41	0.62

下载: 导出CSV

| 显示表格

在坝体施工期结束后,首先模拟了为期55个月的坝体堆石体流变过程,期间库水位在158～182.2 m发生周期性波动。此后,保持库水位158 m,模拟典型夏季1 d的条件,开展了热-力耦合的计算分析。

3.2 面板温度计算结果

图5为计算所得12∶00时面板温度的计算结果。由图5（a）可见,计算面板外表面温度在水位处有明显的分界线,水下面板被约束为水温20℃,水上面板温度较高,表面温度的均值达51.6℃。

图 5 面板的温度分布结果（12∶00时刻）

Figure 5. Spatial distributions of the temperature at the concrete face at the instant of 12∶00

下载: 全尺寸图片幻灯片

图5（b）给出了在河谷中央位置的剖面上,面板温度沿厚度方向的分布。可见,水上面板的温度沿厚度方向由外向内急剧降低,表明太阳辐射所致温升主要集中在面板表层。在图5（b）中局部放大图的右侧,绘制了高程158.9 m处温度沿厚度方向的分布曲线。可见,升温区域主要集中在面板紧邻外表面较浅的薄层中,其中表面1/5厚度内的升温幅度可达到总升温温差的50%以上,表明面板温度随面板厚度方向的幅度衰减十分剧烈。这种分布规律和天生桥一级面板坝的现场观测数据基本一致。这同时也表明,在研究与太阳热辐射相关的问题时,在面板厚度方向划分足够多的单元层数是十分必要的。

图6为计算所得面板最高温度随时间的变化曲线。作为对比,图中同时也给出了太阳热辐射流量 $q_{s}$ 的变化曲线。可见,12∶00为太阳辐射最强的时刻,而面板的最高温度出现在13∶00—14∶00,这是一种典型的相位滞后现象,表明了开展瞬态热传导计算的必要性。

图 6 面板最高温度的时变曲线

Figure 6. Temporal evolution of maximum temperature of concrete face

下载: 全尺寸图片幻灯片

3.3 面板挤压应力计算结果

图7给出了流变结束后尚未考虑温度变化时计算所得面板坝轴向应力的分布。由计算结果可见,由于面板轴向水平位移表现为两岸向河谷中央变形,计算所得面板中间大部分区域属于坝轴向压应力区,左右两岸面板则出现一定范围的拉应力区。面板最大压应力的值为14.9 MPa,发生在河谷中央靠近二期面板顶部的部位。计算结果符合一般情况下高面板坝面板应力的一般规律。需要特别强调的是,在不考虑温度应力的情况下,面板坝轴向应力的最大值并未发生在河谷中央面板的顶部部位。

图 7 不考虑太阳热辐射的面板坝轴向应力计算结果

Figure 7. Distribution of extrusion stress of concrete face under initial temperature

下载: 全尺寸图片幻灯片

图8（a）给出了考虑太阳热辐射13∶00时刻,面板上表面坝轴向应力的计算结果。将其和图7所示的计算结果对比可以发现,水下部分面板的应力基本保持不变。但由于太阳热辐射导致温度升高,使得水上部分面板的轴向应力显著增大。此时,面板坝轴向应力的最大值发生在河谷中部面板的顶部部位,最大值达22.3 MPa,该处应力的最大增幅达到13.3 MPa。需特别指出的是,此时最大坝轴向挤压应力的发生位置,与实际工程中最先发生面板挤压破损的位置相一致。这表明夏季由太阳热辐射所导致的面板温度的升高,应该是使面板发生挤压破损的重要原因之一。

图 8 考虑太阳热辐射的面板坝轴向应力计算结果

Figure 8. Distribution of extrusion stress of concrete face considering solar radiation effect

下载: 全尺寸图片幻灯片

图8（b）给出了计算所得13∶00时刻面板轴向应力沿厚度的分布。由图可见,水上部分面板挤压应力的高值主要分布在面板上表面较薄的一层中,这和图5（b）所示温度的分布规律相一致。这种规律也与实际工程中面板挤压破损主要发生在表层的现象具有高度的相似性。这再次表明,挤压破损应与夏季太阳热辐射所致面板表层温度的升高密切相关。

可将太阳热辐射温度变化导致面板坝轴向压应力增大的原因,概括为两个方面：①轴向挤压作用。面板温度升高导致面板发生热膨胀,使得面板发生坝轴向的挤压作用,从而使得面板坝轴向压应力增大。②温度沿面板厚度变化不均所导致的自应力作用。如前所述,当面板在环境温度和太阳辐射作用下升温时,升温区域主要集中在面板上表面的薄层中,面板下部温度变化不大。在这种情况下,面板表层发生热膨胀,而底部面板则会起到约束的作用,从而在面板上表层区域产生自应力,进一步增大坝轴向的挤压应力。

图9给出了河谷中心断面,不同时刻面板上表面挤压应力沿高程的分布。对照图6温度的计算结果可知,6∶00—13∶00,面板温度处于持续上升期,面板的坝轴向挤压应力也相应持续增大,在13∶00达到最大值22.3 MPa。在出现最大挤压应力的位置处,对应的应力最大增幅为13.3 MPa。13∶00—18∶00,如图9（b）所示,面板温度开始回降,挤压应力也相应逐渐减小。

图 9 不同时刻面板挤压应力沿高程分布

Figure 9. Vertical distribution of extrusion stress of concrete face at different time

下载: 全尺寸图片幻灯片

4. 结论

本文探讨和分析了夏季太阳热辐射温度应力对面板堆石坝混凝土面板坝轴向挤压应力的影响。主要得到5点结论。

（1）采用基于对偶mortar元的计算接触力学方法,推导了可用于非协调网格的非稳定温度场求解格式,发展了考虑温度应力及非线性接触的热-力耦合计算方法,自主开发了相应的有限元数值模拟程序。

（2）基于200 m级理想高面板堆石坝,进行了考虑夏季太阳辐射作用的热力耦合精细化计算分析。为了能准确分析温度及应力沿面板厚度方向的变化,在面板厚度方向划分了6层单元。

（3）计算结果表明,在夏季太阳热辐射作用下,水上面板的上表面可发生较大的温度升高,最大温度可达51.6℃,且升温区域主要集中在面板紧邻上表面的薄层中。计算结果和天生桥一级面板坝现场监测结果的规律总体一致。

（4）太阳热辐射温度升高,会使水上面板的坝轴向挤压应力显著增大,最大值可达22.3 MPa,发生在河谷中部面板的顶部,且挤压应力的高值主要分布在面板上表面的薄层中。这些特点与实际发生面板挤压破损的现象相似,表明夏季太阳热辐射所致的面板温度应力,应是使面板发生挤压破损的重要原因之一。

（5）太阳热辐射导致面板上表面挤压应力增大的作用机制可概括为两个方面：①面板热膨胀导致面板发生轴向挤压；②温度沿面板厚度分布不均所产生的自应力。

图 1 水热侵蚀作用下冻土孔隙扩展示意图

Figure 1. Diagram of pore expansion under fluvial-thermal erosion

下载: 全尺寸图片幻灯片

图 2 双移动边界示意图

Figure 2. Schematic diagram of double-moving boundary

下载: 全尺寸图片幻灯片

图 3 试验系统布置图

Figure 3. Layout of testing apparatus

下载: 全尺寸图片幻灯片

图 4 测温点布置图

Figure 4. Layout of measuring points for temperature

下载: 全尺寸图片幻灯片

图 5 两种情况下测点温度随时间变化的典型特征曲线

Figure 5. Characteristic curves of temperature change with time

下载: 全尺寸图片幻灯片

图 6 典型曲线的融化、潜蚀发生时间点判定示例

Figure 6. Example of determining time point of melting and sub-erosion in typical curves

下载: 全尺寸图片幻灯片

图 7 不同水头差下双边界随时间移动情况

Figure 7. Double boundaries moving with time under different heads

下载: 全尺寸图片幻灯片

图 8 潜蚀初期线性速率v_R0随水头差Δh变化情况

Figure 8. Variation of linear velocity v_R0 with head difference Δh at initial stage of submergence erosion

下载: 全尺寸图片幻灯片

图 9 不同初始温度下双边界随时间变化情况

Figure 9. Variation of double boundaries with time under different initial temperatures

下载: 全尺寸图片幻灯片

图 10 不同冻土含冰量下潜蚀边界(内孔边界)随时间移动情况

Figure 10. Movement of sub-erosion boundary (inner hole boundary) with time under different ice contents of frozen soil

下载: 全尺寸图片幻灯片

表 1 试验分组情况

Table 1 Test groups

组别	水头差Δh/m	初始温度T₀/℃	饱和度S_r
1-1, 2, 3	0.1, 0.15, 0.2	-15	0.92
1-4, 5, 6	0.25, 0.3, 0.35
1-7, 8, 9	0.4, 0.5, 0.6,
1-10, 11	0.8, 1.0
1-12, 13	1.2, 1.4
2-1, 2, 3	0.2	-5, -15, -25	0.92
3-1, 2	0.2	-15	0.8, 0.92
3-3, 4	0.2	-15	1.0, 1.1

下载: 导出CSV

参考文献(9)

[1]	陈瑞杰, 程国栋, 李述训, 等. 人工地层冻结应用研究进展和展望[J]. 岩土工程学报, 2000, 22(1): 40–44. http://cge.nhri.cn/cn/article/id/10448 CHEN Rui-jie, CHENG Guo-dong, LI Shu-xun, et al. Development and prospect of research on application of artificial ground freezing[J]. Chinese Journal of Geotechnical Engineering, 2000, 22(1): 40–44. (in Chinese) http://cge.nhri.cn/cn/article/id/10448
[2]	白云, 汤竞. 软土地下工程的风险管理[J]. 地下空间与工程学报, 2006, 2(1): 21–27. https://www.cnki.com.cn/Article/CJFDTOTAL-BASE200601004.htm BAI Yun, TANG Jing. Risk management for underground project in soft soils[J]. Chinese Journal of Underground Space and Engineering, 2006, 2(1): 21–27. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-BASE200601004.htm
[3]	胡向东, 白楠, 李鸿博. 圣彼得堡地铁1号线区间隧道事故分析[J]. 隧道建设, 2008, 28(4): 418–422. https://www.cnki.com.cn/Article/CJFDTOTAL-JSSD200804007.htm HU Xiang-dong, BAI Nan, LI Hong-bo. Analysis on tunnel accident on line 1 of saint Petersburg metro[J]. Tunnel Construction, 2008, 28(4): 418–422. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JSSD200804007.htm
[4]	HU X D, GUO W, ZHANG L Y, et al. Application of liquid nitrogen freezing to recovery of a collapsed shield tunnel[J]. Journal of Performance of Constructed Facilities, 2014, 28(4): 04014002. doi: 10.1061/(ASCE)CF.1943-5509.0000477
[5]	李柱和. 临江地铁冻结工程事故机理及修复技术研究[D]. 北京: 中国矿业大学(北京), 2013. LI Zhu-he. Research on the Mechanics About Accident of Subway Engineering with Refrigeration Construction Technology and Its Repairing Thechnology[D]. Beijing: China University of Mining & Technology, Beijing, 2013. (in Chinese)
[6]	HU X D, WU Y H, GUO W. Forensic investigation of a water inrush accident during ground freezing recovery work[J]. Journal of Performance of Constructed Facilities, 2019, 33(2): 04019017. doi: 10.1061/(ASCE)CF.1943-5509.0001290
[7]	吴元昊, 胡向东. 柱坐标系下冻土内部热流侵蚀相变问题的求解[J]. 煤炭学报, 2020, 45(2): 660–666. https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB202002015.htm WU Yuan-hao, HU Xiang-dong. Solutions to phase change problem due to internal thermal erosion of frozen soil[J]. Journal of China Coal Society, 2020, 45(2): 660–666. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-MTXB202002015.htm
[8]	李广信, 周晓杰. 土的渗透破坏及其工程问题[J]. 工程勘察, 2004, 32(5): 10–13, 52. https://www.cnki.com.cn/Article/CJFDTOTAL-GCKC200405002.htm LI Guang-xin, ZHOU Xiao-jie. Seepage failure of soil and its problems in engineering[J]. Geotechnical Investigation & Surveying, 2004, 32(5): 10–13, 52. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCKC200405002.htm
[9]	肖朝昀. 人工地层冻结冻土帷幕形成与解冻规律研究[D]. 上海: 同济大学, 2007. XIAO Zhao-yun. Research on the Forming and Thawing of Frozen Soil Walls by Artificial Ground Freezing Method[D]. Shanghai: Tongji University, 2007. (in Chinese)

施引文献(6)

期刊类型引用(4)

1.	姚可夫，田始光，漆一宁，周世华，黄耀英，苏怀智. 高海拔区特征环境驱动下混凝土坝服役性能研究进展. 水利学报. 2023(06): 717-728 . 百度学术
2.	屈春来，吴照正，程方. 太阳辐射下矩形渡槽热力耦合特性分析. 水电能源科学. 2022(07): 159-162+145 . 百度学术
3.	陈生水. 高土石坝变形破坏过程预测理论和防控技术创新. 岩土工程学报. 2022(07): 1211-1219 . 本站查看
4.	汪洋，陈文化. 基于裂隙形状函数的自然环境高温下花岗岩裂隙尖部非线性温度场. 岩土力学. 2022(S1): 267-274 . 百度学术