混凝土面板堆石坝面板太阳热辐射温度应力计算分析

张宗亮; 张丙印; 周墨臻; 尹文杰; 冯业林

doi:10.11779/CJGE202111001

混凝土面板堆石坝面板太阳热辐射温度应力计算分析 English Version

1.
中国电建集团昆明勘测设计研究院有限公司,云南　昆明　650051
2.
清华大学水沙科学与水利水电工程国家重点实验室,北京　100084
3.
北京交通大学城市地下工程教育部重点实验室,北京　100044

基金项目:

国家重点研发计划项目 2017YFC0404802

国家自然科学基金面上项目 51979143

国家自然科学基金青年科学基金项目 51808034

详细信息

作者简介:
张宗亮(1963— ),男,全国工程勘察设计大师,博士,主要从事水利水电工程勘察、设计与科研工作。E-mail：zhang_zl@powerchina.cn

中图分类号: TV641.43
计量
- 文章访问数: 306
- HTML全文浏览量: 35
- PDF下载量: 208
出版历程
- 收稿日期: 2021-03-03
- 网络出版日期: 2022-12-01
- 刊出日期: 2021-10-31

Numerical analysis of thermal stress induced by solar radiation in concrete-faced rockfill dams

1.
Power China Kunming Engineering Corporation Limited, Kunming 650051, China
2.
State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China
3.
Key Laboratory for Urban Underground Engineering of Ministry of Education, Beijing Jiaotong University, Beijing 100044, China

摘要

摘要: 采用基于对偶mortar元的计算接触力学方法,推导了可用于非协调网格的非稳定温度场求解格式,发展了考虑温度应力及非线性接触的热-力耦合计算方法,自主开发了相应的有限元数值模拟程序。对一理想高面板堆石坝进行了考虑夏季太阳辐射作用的热力耦合精细化计算分析。计算结果表明,在夏季太阳热辐射作用下,水上面板的上表面可发生较大的温度升高,最大温度可达51.6℃,并造成面板坝轴向挤压应力显著增大,最大可达22.3 MPa。此时,面板最大挤压应力发生在河谷中部面板的顶部,且挤压应力的高值主要分布在面板上表面的薄层中。这些特点与实际工程发生面板挤压破损的现象相似,表明夏季太阳热辐射所致的面板温度应力,是使面板发生挤压破损的重要原因之一。
- 堆石坝 /
- 面板挤压破损 /
- 太阳热辐射 /
- 对偶mortar有限元 /
- 温度应力
Abstract: Based on the dual mortar method within the framework of computational contact mechanics, the transient heat conduction theory with nonconforming meshes is derived. A new thermo-mechanical coupled method is proposed to account for the temperature-induced stress and nonlinear contact behavior in concrete-faced rockfill dams (CFRDs). It is implemented numerically in the authors' in-house finite element code and then used to simulate an ideal CFRD under the solar radiation effect. The numerical results indicate that the solar radiation in summer can cause a significant temperature increase to the concrete face above the reservoir water. The maximum temperature reaches up to 51.6 ℃ and hence contributes a considerable additional value to the axial extrusion stress of the concrete face. The maximum stress is 22.3 MPa and occurs at the central valley near the dam top. The high values of extrusion stress are located at a thin layer near the outer surface. The numerical phenomenon coincides with the common features of extrusion damage from the practical CFRDs, confirming that the thermal stress induced by the summer solar radiation is one of the major factors causing extrusion damage.
- concrete-faced rockfill dam /
- extrusion damage /
- solar radiation /
- dual mortar finite element method /
- thermal stress

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. Extrusion damages and calculated extrusion stresses of concrete face in Tianshengqiao-1 CFRD

下载: 全尺寸图片幻灯片

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

Figure 2. Three boundary conditions in thermal analysis

下载: 全尺寸图片幻灯片

图 3 理想坝的材料分区

Figure 3. Material zones and filling stages of ideal dam

下载: 全尺寸图片幻灯片

图 4 三维计算网格

Figure 4. Three-dimensional computational mesh

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

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

Figure 6. Temporal evolution of maximum temperature of concrete face

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

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

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

下载: 全尺寸图片幻灯片

表 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}$ $α_{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

表 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

参考文献(25)

[1]	曹克明, 徐建军. 超高面板坝面板临界挠度探讨和设计改进[J]. 水力发电, 2008, 34(11): 98-102. doi: 10.3969/j.issn.0559-9342.2008.11.028 CAO Ke-ming, XU Jian-jun. Discussions on critical deflection of face slab and its design improvement for super-high CFRD[J]. Water Power, 2008, 34(11): 98-102. (in Chinese) doi: 10.3969/j.issn.0559-9342.2008.11.028
[2]	郦能惠. 高混凝土面板堆石坝设计理念探讨[J]. 岩土工程学报, 2007, 29(8): 1143-1150. doi: 10.3321/j.issn:1000-4548.2007.08.005 LI Neng-hui. New concept of design for high concrete face rockfill dams[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(8): 1143-1150. (in Chinese) doi: 10.3321/j.issn:1000-4548.2007.08.005
[3]	马洪琪, 曹克明. 超高面板坝的关键技术问题[J]. 中国工程科学, 2007, 9(11): 4-10. doi: 10.3969/j.issn.1009-1742.2007.11.001 MA Hong-qi, CAO Ke-ming. Key technical issues related to super-high concrete slab dam[J]. Engineering Sciences, 2007, 9(11): 4-10. (in Chinese) doi: 10.3969/j.issn.1009-1742.2007.11.001
[4]	罗亮, 陈烨, 仲洪滔. 土工膜在大坝面板挤压破坏临时处理工程中的应用[J]. 大坝与安全, 2013(2): 48-51. doi: 10.3969/j.issn.1671-1092.2013.02.012 LUO Liang, CHEN Ye, ZHONG Hong-tao. Application of geomembrane in temporary treatment of extrusion damage of dam face slab[J]. Dam and Safety, 2013(2): 48-51. (in Chinese) doi: 10.3969/j.issn.1671-1092.2013.02.012
[5]	郦能惠, 王君利, 米占宽, 等. 高混凝土面板堆石坝变形安全内涵及其工程应用[J]. 岩土工程学报, 2012, 34(2): 193-201. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201202003.htm LI Neng-hui, WANG Jun-li, MI Zhan-kuan, et al. Connotation of deformation safety of high concrete face rockfill dams and its application[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(2): 193-201. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201202003.htm
[6]	徐泽平, 郭晨. 高面板堆石坝面板挤压破坏问题研究[J]. 水力发电, 2007, 33(9): 80-84. doi: 10.3969/j.issn.0559-9342.2007.09.029 XU Ze-ping, GUO Chen. Research on the concrete face slab rupture of high CFRD[J]. Water Power, 2007, 33(9): 80-84. (in Chinese) doi: 10.3969/j.issn.0559-9342.2007.09.029
[7]	邹德高, 尤华芳, 孔宪京, 等. 接缝简化模型及参数对面板堆石坝面板应力及接缝位移的影响研究[J]. 岩石力学与工程学报, 2009, 28(增刊1): 3257-3263. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2009S1103.htm ZOU De-gao, YOU Hua-fang, KONG Xian-jing, et al. Research on joint simplified model and effects of joint parameters on panel stress and joint displacements of faced rockfill dam[J]. Chinese Journal of Rock Mechanics and Engineering, 2009, 28(S1): 3257-3263. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2009S1103.htm
[8]	周墨臻, 张丙印, 张宗亮, 等. 超高面板堆石坝面板挤压破坏机理及数值模拟方法研究[J]. 岩土工程学报, 2015, 37(8): 1426-1432. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201508014.htm ZHOU Mo-zhen, ZHANG Bing-yin, ZHANG Zong-liang, et al. Mechanisms and simulation methods for extrusion damage of concrete faces of high concrete-faced rockfill dams[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(8): 1426-1432. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201508014.htm
[9]	刘光廷, 胡昱, 焦修刚, 等. 高面板堆石坝面板应力规律分析及改善应力状态的对策[J]. 水利学报, 2006, 37(2): 135-140. doi: 10.3321/j.issn:0559-9350.2006.02.002 LIU Guang-ting, HU Yu, JIAO Xiu-gang, et al. Stress analysis of concrete slab in high CFRD and measures for improving the stress state[J]. Journal of Hydraulic Engineering, 2006, 37(2): 135-140. (in Chinese) doi: 10.3321/j.issn:0559-9350.2006.02.002
[10]	王瑞骏, 王党在, 陈尧隆. 寒潮冷击作用下堆石坝混凝土面板温度应力研究[J]. 水力发电学报, 2004, 23(6): 45-49. doi: 10.3969/j.issn.1003-1243.2004.06.010 WANG Rui-jun, WANG Dang-zai, CHEN Yao-long. Research on thermal stress of concrete slabs of rockfill dam under condition of cold wave attack[J]. Journal of Hydroelectric Engineering, 2004, 23(6): 45-49. (in Chinese) doi: 10.3969/j.issn.1003-1243.2004.06.010
[11]	张国新, 张丙印, 王光纶. 混凝土面板堆石坝温度应力研究[J]. 水利水电技术, 2001, 32(7): 1-5, 62. https://www.cnki.com.cn/Article/CJFDTOTAL-SJWJ200107000.htm ZHANG Guo-xin, ZHANG Bing-yin, WANG Guang-lun. Study of thermal stress in the concrete slab of a concrete face rock-fill dam[J]. Water Resources and Hydropower Engineering, 2001, 32(7): 1-5, 62. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SJWJ200107000.htm
[12]	程嵩, 张嘎, 张建民, 等. 有挤压墙面板堆石坝的面板温度应力分析及改善措施研究[J]. 工程力学, 2011, 28(4): 76-81. https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201104014.htm CHENG Song, ZHANG Ga, ZHANG Jian-min, et al. The analysis and amendatory measures of temperature stress of face slab of a concrete-faced rockfill dam with extrusion wall[J]. Engineering Mechanics, 2011, 28(4): 76-81. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201104014.htm
[13]	颉志强, 张振杰, 董芸, 等. 基于数值仿真的堆石坝面板施工期温控防裂方法研究[J]. 水利水电技术, 2018, 49(增刊1): 33-41. https://www.cnki.com.cn/Article/CJFDTOTAL-SJWJ2018S1006.htm XIE Zhi-qiang, ZHANG Zhen-jie, DONG Yun, et al. Numerical simulation-based study of method for temperature control and crack prevention of concrete face slab of rock-fill dam during construction[J]. Water Resources and Hydropower Engineering, 2018, 49(S1): 33-41. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SJWJ2018S1006.htm
[14]	王子健, 刘斯宏, 李玲君, 等. 公伯峡面板堆石坝面板裂缝成因数值分析[J]. 水利学报, 2014, 45(3): 343-350. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201403012.htm WANG Zi-jian, LIU Si-hong, LI Ling-jun, et al. Numerical analysis of the causes of slab’s cracks on Gongboxia face rockfill dam[J]. Journal of Hydraulic Engineering, 2014, 45(3): 343-350. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201403012.htm
[15]	宋文晶, 高莲士, 吕明治, 等. 张河湾水电站上水库沥青混凝土面板坝应力变形分析[J]. 水力发电学报, 2007, 26(4): 82-85. https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB200704016.htm SONG Wen-jing, GAO Lian-shi, LÜ Ming-zhi, et al. Stress and deformation FEM analysis of Zhanghewan asphalt CFRD[J]. Journal of Hydroelectric Engineering, 2007, 26(4): 82-85. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB200704016.htm
[16]	石北啸, 蔡正银, 陈生水. 温度变化对堆石料变形影响的试验研究[J]. 岩土工程学报, 2016, 38(增刊2): 299-305. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2016S2049.htm SHI Bei-xiao, CAI Zheng-yin, CHEN Sheng-shui. Experiments on influence of temperature on deformation of rock fills[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(S2): 299-305. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2016S2049.htm
[17]	徐泽平. 混凝土面板堆石坝关键技术与研究进展[J]. 水利学报, 2019, 50(1): 62-74. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201901008.htm XU Ze-ping. Research progresses and key technologies of CFRD construction[J]. Journal of Hydraulic Engineering, 2019, 50(1): 62-74. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201901008.htm
[18]	PUSO M A, LAURSEN T A. A mortar segment-to-segment frictional contact method for large deformations[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193(45/46/47): 4891-4913.
[19]	WOHLMUTH B I. A mortar finite element method using dual spaces for the Lagrange multiplier[J]. SIAM Journal on Numerical Analysis, 2000, 38(3): 989-1012.
[20]	POPP A, SEITZ A, GEE M W, et al. Improved robustness and consistency of 3D contact algorithms based on a dual mortar approach[J]. Computer Methods in Applied Mechanics and Engineering, 2013, 264: 67-80.
[21]	FENG S Z, CUI X Y, LI G Y. Transient thermal mechanical analyses using a face-based smoothed finite element method (FS-FEM)[J]. International Journal of Thermal Sciences, 2013, 74: 95-103.
[22]	朱伯芳. 大体积混凝土施工过程中受到的日照影响[J]. 水力发电学报, 1999, 18(3): 35-41. https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB903.004.htm ZHU Bo-fang. Influence of solar radiation on temperature of mass concrete in the process of construction[J]. Journal of Hydroelectric Engineering, 1999, 18(3): 35-41. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SFXB903.004.htm
[23]	ZHOU Mo-zhen, ZHANG Bing-yin, PENG Chong, et al. Three-dimensional numerical analysis of concrete-faced rockfill dam using dual-mortar finite element method with mixed tangential contact constraints[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2016, 40(15): 2100-2122.
[24]	KUNISCH K, RÖSCH A. Primal-dual active set strategy for a general class of constrained optimal control problems[J]. SIAM Journal on Optimization, 2002, 13(2): 321-334.
[25]	李国英, 米占宽, 傅华, 等. 混凝土面板堆石坝堆石料流变特性试验研究[J]. 岩土力学, 2004, 25(11): 1712-1716. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200411008.htm LI Guo-ying, MI Zhan-kuan, FU Hua, et al. Experimental studies on rheological behaviors for rockfills in concrete faced rockfill dam[J]. Rock and Soil Mechanics, 2004, 25(11): 1712-1716. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200411008.htm

施引文献(6)

期刊类型引用(4)

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