面板堆石坝混凝土面板配筋计算方法研究

魏匡民; 周恒; 米占宽; 李国英; 邓曌

doi:10.11779/CJGE20220757

面板堆石坝混凝土面板配筋计算方法研究 English Version

魏匡民^{1, 3,},
周恒²,
米占宽^{1, 3},
李国英^{1, 3},
邓曌¹

1.
南京水利科学研究院, 江苏南京 210024
2.
中国电建集团西北勘测设计研究院有限公司, 陕西西安 710065
3.
水利部水库大坝安全重点实验室, 江苏南京 210029

基金项目:

国家重点研发计划 2021YFC3090100

中央级公益性科研院所基本科研业务费项目 Y322003

中央级公益性科研院所基本科研业务费项目 Y320005

国家自然科学基金项目 52109094

江苏省水利科技项目 2018019

详细信息

作者简介:
魏匡民(1985—)，男，工学博士，高级工程师，主要从事土石坝安全保障理论与技术方面的研究工作。E-mail: kmwei@nhri.cn

中图分类号: TU43;TV641
计量
- 文章访问数: 452
- HTML全文浏览量: 78
- PDF下载量: 91
出版历程
- 收稿日期: 2022-06-13
- 网络出版日期: 2023-02-19
- 刊出日期: 2023-05-31

Reinforcement calculation of concrete slab in concrete faced rockfill dams

WEI Kuangmin^{1, 3,},
ZHOU Heng²,
MI Zhankuan^{1, 3},
LI Guoying^{1, 3},
DENG Zhao¹

1.
Nanjing Hydraulic Research Institute, Nanjing 210024, China
2.
Northwest Engineering Corporation Limited, Xi'an 710065, Chian
3.
Key Laboratory of Reservoir and Dam Safety of the Ministry of Water Resources, Nanjing 210029, China

摘要

摘要: 混凝土面板应力变形安全是面板堆石坝工程的关键，然而传统面板坝分析方法大多未考虑面板的配筋效应，实际面板配筋方案制定缺乏理论依据。提出了基于限裂理念的混凝土面板配筋计算方法：采用非协调网格计算技术精细模拟混凝土面板应力，采用钢筋埋置单元模拟钢筋加固作用，采用损伤模型模拟混凝土材料损伤开裂行为，通过引入与几何尺寸、钢筋应力相关的裂缝宽度公式实现混凝土面板限裂配筋计算。提出的方法可定量研究面板限裂配筋，为复杂条件下面板坝配筋设计提供理论依据。文中以正在施工建设的大石峡特高面板坝为例，初步研究了该坝的面板配筋方案。
- 面板堆石坝 /
- 有限元 /
- 面板配筋计算 /
- 精细模拟 /
- 损伤模型
Abstract: The stress and deformation safety of the concrete slab is crucial to the concrete faced rockfill dams (CFRDs), however, most of the traditional methods do not consider the effects of the reinforcement in slab, therefore, the calculated results can not provide enough information to the design of slab reinforcement. In this study, a method for calculating the slab reinforcement based on the concept of crack width limit is proposed: the non-coordinated grid computing technology is used to refinedly simulate the stress of the concrete slab, the embedded steel element is used to simulate the reinforcement effects, and the damage model is used to simulate the cracking of the slab concrete. By introducing a formula for crack width related to geometric size and rebar stress, and setting the crack width limit, the reinforcement calculation is realized. The proposed method can be used to quantitatively study the slab reinforcement and provide a theoretical basis for reinforcement design under complex conditions. The Dashixia CFRD under construction is taken as an example and its slab reinforcement design is preliminarily studied.
- concrete faced rockfill dam /
- finite element /
- reinforcement calculation /
- refined simulation /
- damage model

HTML全文

0. 引言

在基坑工程中，地下水回灌常被用于控制或消除基坑降水对周围环境的影响。地下水回灌可以分为定流量回灌和定水头回灌，其中定水头回灌常用于低渗透性含水层的水位抬升^[1-3]。最早的定水头井流模型是由Jacob等^[4]建立的。该模型描述了水平无限延伸的均质承压含水层中完整井定降深条件下的径向水流运动。在此基础上，国内外众多学者考虑含水层越流、表皮效应、非完整井、含水层类型等因素的影响，对Jacob等模型进行了改进并推导了相应的解析解或半解析解^[1-3]。由于地下水回灌通常被视为地下水开采的逆过程，上述单井定降深井流模型可同样为定水头回灌条件下的含水层渗流响应分析提供参考。然而，由于相反的水流流动方向改变了井周的物理条件，不同于抽汲地下水在井周形成的渗透增高地带，回灌地下水携带的细颗粒会随着水流速度下降逐渐沉积而在井周形成一个渗透衰减地带^[5]。与此同时，回灌水中溶解的气泡、回灌水与地下水或含水层介质间的化学反应以及滤管周围微生物的生长繁殖均可能进一步加剧井周含水层渗透性衰减。上述物理、生物与化学过程导致的渗透衰减使得基于恒定渗透系数的抽水井流模型不再适用于分析回灌井流问题。

为考虑回灌过程中渗透衰减效应，何满潮等^[6]通过对地热水回灌的测试数据分析，建立了考虑井周物理与化学阻塞效应的渗透系数衰减模型。随后，何满潮等^[5]将该服从指数衰减的渗透模型引入经典的Theis模型^[7]，推导了考虑渗透系数动态衰减的单井回灌公式，并分析了渗透衰减对于回灌井周渗流场的影响。基于上述单井回灌公式，闫峭等^[8]利用改进遗传算法研究了考虑渗透系数衰减的水文地质参数反演问题。上述研究均只考虑了定流量回灌，并未涉及定水头回灌。对于由回灌压力控制的回灌井流问题，李炯等^[9]将何满潮等^[6]提出的渗透衰减方程引入至Jacob等^[4]模型，并利用Weber变换得到了模型的解析解，但该模型未考虑回灌井非完整性的影响。此外，上述研究均假设整个含水层的渗透系数存在同步衰减，而现实中回灌堵塞通常仅存在于回灌井邻近的柱状包围带^[5]，堵塞对主含水层渗透性的影响非常轻微。

为分析定水头非完整井回灌过程中井周堵塞效应对承压含水层渗流场的影响，本文选用指数形式的渗透方程代替恒定渗透系数反映堵塞作用下的井周含水层渗透衰减行为，建立了考虑井周堵塞效应的非完整井定水头回灌井流模型。通过变量代换、Laplace变换和有限余弦Fourier变换及其逆变换，推导了井流模型在Laplace空间的解，并利用Stehfest数值逆变换方法得到了实时空间内的半解析解，并进一步分析了回灌堵塞诱发的渗透衰减对定水头非完整井回灌渗流场的影响规律。研究结果可以为定水头回灌的工程设计、参数反演、堵塞治理提供理论参考。

1. 问题的提出与求解

1.1 渗透衰减模型

由于物理、生物以及化学因素引起的堵塞在回灌井周围并不是均匀分布的，距离回灌井越近，地层堵塞越严重，渗透系数下降越显著，使得回灌堵塞影响下的井周含水层渗透系数兼有时间和空间变异性。然而，当前试验条件下准确获得径向流动过程中的堵塞时空发展规律及其导致的渗透时空变异规律是非常困难的。同时，考虑时空变异的渗透衰减方程对于井流模型的解析求解亦提出了更高的要求。为便于求解与分析，本文仅考虑渗透系数的时变效应。需要指出的是，这里的渗透系数表征的是整个井周堵塞区而非某一给定位置的水力传导性能，是对堵塞区实际渗透系数在空间上的平均唯象近似。

基于上述分析，为考虑回灌堵塞对含水层导水特性的影响，本文假设井周含水层在水平和竖直方向上的渗透系数均服从Bianchi等^[10]基于室内水平砂柱试验建立的反映堵塞作用下渗透衰减问题的概括性方程。该方程具有如下形式：

${K_{\text{r}}}(t) = {K_{{\text{r}},\infty }} + ({K_{{\text{r}},0}} - {K_{{\text{r}},\infty }})\exp ( - {\lambda _0} \cdot t) \;\;\text{，}$

(1a)

${K_{\text{z}}}(t) = {K_{{\text{z}},\infty }} + ({K_{{\text{z}},0}} - {K_{{\text{z}},\infty }})\exp ( - {\lambda _0} \cdot t) \;\;。$

(1b)

式中：K_r和K_z分别为含水层水平和竖直方向上的渗透系数；“0”和“∞”分别表示渗透系数的初始与渐近值； ${\lambda _0}$ 为渗透系数随时间的衰减指数。当K_r = K_z且K_∞ = 0时，上述渗透衰减模型可退化至何满潮等^[6]提出的渗透衰减模型。此外，为便于数学分析，本文假设含水层水平和竖向渗透系数具有相同的衰减速率，即相同的 ${\lambda _0}$ ，且两者之比的初始值和稳定值保持不变，即K_{r, 0}/K_{z, 0} = K_{r, ∞}/K_{z, ∞} = ${\alpha ^2}$ 。此时，在整个渗透衰减过程中K_r(t)与K_z(t)的比值亦为定值 ${\alpha ^2}$ 。

1.2 数学模型

本文选用的非完整井定水头回灌模型位于图 1所示的含水层系统中。该含水层系统由一个均质等厚、侧向无限延伸的水平各向同性承压含水层及其上覆和下伏隔水层组成，是一种工程实践中常见的含水层系统，便于开展相关的野外现场测试和室内模型试验，从而有利于验证理论模型的适用性和准确性。同时，该含水层系统也被众多经典井流模型选用，如Yang等^[3]建立的考虑表皮效应的两区井流模型。鉴于此，本文采用了与Yang等模型几乎完全相同的含水层假设。但不同于Yang等假定井周表皮层渗透系数随时间恒定，本文模型须考虑回灌堵塞造成的井周含水层渗透衰减效应。对于井周含水层回灌堵塞区，参考Yang等^[3]对表皮层的处理方式，本文假设井周堵塞区的径向延伸r_c沿深度保持不变，呈圆柱状分布，如图 1所示。在该柱状包围带内，含水层具有良好的均质性，不同位置的渗透系数受堵塞效应影响而随时间同步衰减。

图 1 承压含水层中非完整井定水头回灌示意

Figure 1. Schematic graph of constant-head injection in confined aquifers using a partially penetrating well

下载: 全尺寸图片幻灯片

基于上述假设，以回灌井轴线与含水层底板交点为原点、回灌井轴线竖直向上为z轴、含水层底板水平径向延伸为r轴建立柱坐标系。分别用式（1a）和式（1b）代替Yang等^[3]模型中水平与竖直方向上的常渗透系数，则考虑堵塞效应的非完整井定水头回灌模型的控制方程可以描述为

${K_{\text{r}}}(t)\left[ {\frac{{{\partial ^2}s}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial s}}{{\partial r}}} \right] + {K_{\text{z}}}(t)\frac{{{\partial ^2}s}}{{\partial {z^2}}} = \frac{S}{B}\frac{{\partial s}}{{\partial t}}\;\;。$

(2)

式中：s为含水层的水头抬升；r和z分别为到回灌井轴线的径向距离和到含水层底板的垂向距离；t为回灌时间；S为含水层储水系数；B为含水层厚度。

回灌开始前，井周含水层初始水头变化为零：

$s(r,z,t) = 0,{\text{ }}t = 0\;\;。$

(3)

回灌过程中，井筒内的平均压力水头恒为s_w：

$\frac{1}{{l - d}}{\left. {\int_d^l {s(r,z,t)} {\text{d}}z} \right|_{r = {r_{\text{w}}}}} = {s_{\text{w}}}{\text{ (}}d \leqslant z \leqslant l) \;\;。$

(4)

式中：r_w为回灌井半径；l和d分别为回灌井过滤器顶端与底端的竖向坐标。

在回灌井内，假设定水头作用下回灌流量Q(t)沿滤管均匀分布，则有

${\left. { - 2{\text{π }}r(l - d){K_{\text{r}}}(t){{\partial s(r,z,t)} \mathord{\left/ {\vphantom {{\partial s(r,z,t)} {\partial r}}} \right. } {\partial r}}} \right|_{r = {r_{\text{w}}}}}{\text{ = }}Q(t){\text{ (}}d \leqslant z \leqslant l) \;\;。$

(5)

回灌过程中，沿井壁套筒流量为零，含水层与下伏和上覆弱透水层的交界面完全隔水，则有

${\left. {{{\partial s(r,z,t)} \mathord{\left/ {\vphantom {{\partial s(r,z,t)} {\partial r}}} \right. } {\partial r}}} \right|_{r = {r_{\text{w}}}}} = {\text{0 (}}z \leqslant d,{\text{ }}z \geqslant l) \;\;\text{，}$

(6)

${{\partial s(r,z,t)} \mathord{\left/ {\vphantom {{\partial s(r,z,t)} {\partial z}}} \right. } {\partial z}} = 0\;\;\;{\text{ (}}z = 0)\;\;\text{，}$

(7)

${{\partial s(r,z,t)} \mathord{\left/ {\vphantom {{\partial s(r,z,t)} {\partial z}}} \right. } {\partial z}} = 0\;\;\;{\text{ (}}z = B)\;\;。$

(8)

回灌过程中，井周堵塞区与主含水层的交界面上的水头抬升如下式所示：

$s\left(r_{\mathrm{c}}, z, t\right)=s_{\mathrm{c}}(z, t) \quad(0 \leqslant z \leqslant B) \;\;。$

(9)

式中，s_c(z, t)为回灌流量等于Q(t)且不考虑井周堵塞效应时求得的r = r_c处的水头抬升，其表达式与Feng等^[1]推导的变流量条件下承压含水层中水头变化公式相同。式（9）实际上忽略了井周堵塞对主含水层的渗流场的影响。根据Yeh等^[11]对考虑表皮效应的两区井流模型的分析结果可知，试井周围的低渗透区域仅会在试验初期对主含水层渗流场有轻微影响。随着井周渗流场在很短时间内达到拟稳态，该影响仅会在试验初期带来非常微小的误差，是可以接受的。

1.3 问题的求解

式（2）~（9）构成了考虑井周堵塞效应的非完整井定水头回灌井流模型的定解问题。为便于求解，首先对定解问题作如下变量代换：

$r^{\prime} = {r \mathord{\left/ {\vphantom {r \alpha }} \right. } \alpha } \;\;\text{，}$

(10)

$t^{\prime} = {K_{{\text{z}},\infty }}t/{K_{{\text{z}},0}} + ({K_{{\text{z}},0}} - {K_{{\text{z}},\infty }})(1 - {{\text{e}}^{ - {\lambda _0}t}})/({\lambda _0}{K_{{\text{z}},0}}) \;\;\text{，}$

(11)

式中， $\alpha = \sqrt {{{{K_{{\text{r,0}}}}} \mathord{\left/ {\vphantom {{{K_{{\text{r,0}}}}} {{K_{{\text{z,0}}}}}}} \right. } {{K_{{\text{z,0}}}}}}} = \sqrt {{{{K_{{\text{r,}}\infty }}} \mathord{\left/ {\vphantom {{{K_{{\text{r,}}\infty }}} {{K_{{\text{z,}}\infty }}}}} \right. } {{K_{{\text{z,}}\infty }}}}}$ 。经过坐标变换之后，含水层的水头抬升s(r, z, t)用s^'(r^', z, t^')表示，则定解问题转换为

$\frac{{{\partial ^2}s^{\prime}}}{{\partial {{r^{\prime}}^2}}} + \frac{1}{{r^{\prime}}}\frac{{\partial s^{\prime}}}{{\partial r^{\prime}}} + \frac{{{\partial ^2}s^{\prime}}}{{\partial {z^2}}} = \frac{S}{{B{K_{{\text{z}},0}}}}\frac{{\partial s^{\prime}}}{{\partial t^{\prime}}} \;\;\text{，}$

(12)

$s^{\prime}\left( {r^{\prime},z,t^{\prime}} \right) = 0,{\text{ }}t^{\prime} = 0\;\;\text{，}$

(13)

$\frac{1}{{l - d}}{\left. {\int_d^l {s^{\prime}(r^{\prime},z,t^{\prime})} {\text{d}}z} \right|_{r^{\prime} = {{r}_{\text{w}}^{\prime}}}} = {s_{\text{w}}}\;\;\;{\text{ (}}d \leqslant z \leqslant l) \;\;\text{，}$

(14)

${\left. {\frac{{\partial s^{\prime}(r^{\prime},z,t^{\prime})}}{{\partial r^{\prime}}}} \right|_{r^{\prime} = {{r}_{\text{w}}^{\prime}}}}{\text{ = }} - \frac{{{{Q^{\prime}}_1}(t^{\prime})}}{{2{\text{π }}{{r}_{\text{w}}^{\prime}}(l - d){K_{{\text{r,0}}}}}}\;\;\;{\text{ (}}d \leqslant z \leqslant l) \;\;\text{，}$

(15)

${\left. {{{\partial s^{\prime}(r^{\prime},z,t^{\prime})} \mathord{\left/ {\vphantom {{\partial s^{\prime}(r^{\prime},z,t^{\prime})} {\partial r^{\prime}}}} \right. } {\partial r^{\prime}}}} \right|_{r^{\prime} = {{r}_{\text{w}}^{\prime}}}} = 0\;\;\;{\text{ (}}z \leqslant d,{\text{ }}z \geqslant l) \;\;\text{，}$

(16)

${{\partial s^{\prime}(r^{\prime},z,t^{\prime})} \mathord{\left/ {\vphantom {{\partial s^{\prime}(r^{\prime},z,t^{\prime})} {\partial z}}} \right. } {\partial z}} = 0\;\;\;{\text{ (}}z = 0)\;\;\text{，}$

(17)

${{\partial s^{\prime}{\text{(}}r^{\prime},z,t^{\prime}{\text{)}}} \mathord{\left/ {\vphantom {{\partial s^{\prime}{\text{(}}r^{\prime},z,t^{\prime}{\text{)}}} {\partial z}}} \right. } {\partial z}} = 0\;\;\;{\text{ (}}z = B)\;\;\text{，}$

(18)

$s^{\prime}({r^{\prime}_{\text{c}}},z,t^{\prime}){\text{ = }}{s^{\prime}_{\text{c}}}(z,t^{\prime})\;\;\;{\text{ (}}0 \leqslant z \leqslant B) \;\;。$

(19)

其中，式（15）中 ${Q^{\prime}_1}\left( {t^{\prime}} \right){\text{ = }}{{Q^{\prime}\left( {t^{\prime}} \right){K_{{\text{r,0}}}}} \mathord{\left/ {\vphantom {{Q^{\prime}\left( {t^{\prime}} \right){K_{{\text{r,0}}}}} {{{K^{\prime}}_{\text{r}}}\left( {t^{\prime}} \right)}}} \right. } {{{K^{\prime}}_{\text{r}}}\left( {t^{\prime}} \right)}}$ ，Q^'(t^')和K_r^'(t^')以及式（19）中的 ${s^{\prime}_{\text{c}}}(z,t^{\prime})$ 分别为回灌流量Q(t)、渗透系数方程K_r(t)以及堵塞区外部水头抬升 ${s_{\text{c}}}(z,t)$ 作关于时间的变量代换之后得到的表达式。由于根据式（11）获得t^'关于t的表达式十分困难，使得上述表达式的具体形式同样难以确定。根据 ${s^{\prime}_{\text{c}}}(z,t^{\prime}) = {s_{\text{c}}}(z,t(t^{\prime}))$ ， ${s^{\prime}_{\text{c}}}(z,t^{\prime})$ 可表示为 ${s_{\text{c}}}(z,t^{\prime})$ 与 $\Delta {s_{\text{c}}}(z,t^{\prime})$ 之和，其中 $\Delta {s_{\text{c}}}(z,t^{\prime})$ 的表达式通过非线性拟合求得。考虑到后续求解需要对 $\Delta {s_{\text{c}}}(z,t^{\prime})$ 分别进行关于z和t^'的积分变换，因此 $\Delta {s_{\text{c}}}(z,t^{\prime})$ 的拟合公式宜采用如下可分离变量的形式：

$\Delta {s_{\text{c}}}(z,t^{\prime}) = \Delta {s_{\text{c}}}(z) \cdot \Delta {s_{\text{c}}}(t^{\prime}) \;\;。$

(20a)

同时，根据Laplace变换和有限Fourier变换的定义，当原函数为指数函数时，可以直接求得积分变换后的象函数。因此， $\Delta {s_{\text{c}}}(z)$ 和 $\Delta {s_{\text{c}}}(z)$ 宜采用如下形式的分段指数函数进行拟合：

$\Delta {s_{\text{c}}}(t^{\prime}) = \left[ {U(t^{\prime} - {{t^{\prime}}_{i + 1}}) - U(t^{\prime} - {{t^{\prime}}_i})} \right]\left[ {{A_i} + {A_{i + 1}}\exp ({\alpha _i}t^{\prime})} \right] \;\;\text{，}$

(20b)

$\Delta {s_{\text{c}}}(z) = \left[ {U(z - {z_{j + 1}}) - U(z - {z_j})} \right]\left[ {{C_j} + {C_{j + 1}}\exp ({\varepsilon _j}z)} \right] \;\;。$

(20c)

式中：A_i， ${\alpha _i}$ ，C_j及ε_j为根据 ${s_{\text{c}}}\left[ {z,t(t^{\prime})} \right]$ 与 ${s_{\text{c}}}(z,t^{\prime})$ 的差值序列确定的拟合参数；U(·)为Heaviside跃阶函数，当t < ${t_i}$ 时，U $(t - {t_i})$ = 0，当t ≥ ${t_i}$ 时，U $(t - {t_i})$ =1， ${t^{\prime}_i}$ 和 ${z_j}$ 为分段点。

对定解问题作关于时间 $t^{\prime}$ 的Laplace变换和关于z的有限Fourier余弦变换可得

$\frac{\partial^2 \hat{\bar{s}^{\prime}}}{\partial r^{\prime 2}}+\frac{1}{r^{\prime}} \frac{\partial \hat{\bar{s}^{\prime}}}{\partial r^{\prime}}-\xi(n, p)^2 \hat{\bar{s}^{\prime}}=0\;\;\;,$

(21)

$\frac{1}{{l - d}}\int_d^l {\bar s^{\prime}({{r}_{\text{w}}^{\prime}},z,p)} {\text{d}}z = \frac{{{s_{\text{w}}}}}{p}{\text{ (}}d \leqslant z \leqslant l)\;\;\text{，}$

(22)

$\left.K_{\mathrm{r}, 0} \frac{\partial \hat{\bar{s}^{\prime}}\left(r^{\prime}, n, p\right)}{\partial r^{\prime}}\right|_{r^{\prime}=r_{\mathrm{w}}^{\prime}}=\frac{-\bar{Q}_1^{\prime}(p)}{2 \pi r_{\mathrm{w}}^{\prime}(l-d)} \cdot A(n)\;\;\text{，}$

(23)

$\hat{\bar{s}^{\prime}}({r^{\prime}_{\text{c}}},n,p){\text{ = }}{\hat{\bar{s}^{\prime}}_{\text{c}}}(n,p) + \Delta {\hat{\bar{s}^{\prime}}_{\text{c}}}(n,p) \;\;。$

(24)

式中： $\hat{\bar{s}^{\prime}}(r^{\prime},n,p)$ 和 $\Delta {\hat{\bar{s}^{\prime}}_{\text{c}}}(n,p)$ 分别为水头抬升 $s(r^{\prime},n,p)$ 和水头差 $\Delta {s_{\text{c}}}(n,p)$ 在Laplace-Fourier空间的表达式，符号“-”和“^”分别表示Laplace和Hankel空间的项，p和n分别为Laplace变量和Fourier变量， $\xi {({r}^{\prime },n,p)}^{2}={w}_{n}^{2}+Sp/(B{K}_{\text{z},0})，A(n)=[\mathrm{sin}({w}_{n}l)-$ $\sin ({w_n}d)]/{w_n}$ ， ${w_n} = n{\text{π /}}B$ 。式（21）为零阶修正Bessel方程，通解为

$\begin{array}{l} \hat{\bar{s}}(r^{\prime},n,p) = {C_1}(n,p){{\text{I}}_0}\left[ {{\xi _0}(n,p)r^{\prime}} \right] + {C_2}(n,p){{\text{K}}_0} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {{\xi _0}(n,p)r^{\prime}} \right] \;\;。 \end{array}$

(25a)

式中：待定系数C₁(n, p)和C₂(n, p)根据边界条件式（23）和（24）确定表达式如下所示：

${C}_{1}(n,p)=\frac{\widehat{\overline{{s}}^{\prime}}({{r}}_{\text{c}}^{\prime },n,p){\text{K}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]-\frac{{\overline{{Q}}}_{1}^{\prime}(p)A(n){\text{K}}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]}{2\text{π}{{r}}_{\text{w}}^{\prime }(l-d){K}_{\text{r,0}}{\xi }_{0}(n,p)}}{{\text{I}}_{\text{1}}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{K}}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]+{\text{K}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{I}}_{0}\left[{\xi }_{0}(n,p){{r}^{\prime }}_{}\right]}\text{，}$

(25b)

${C}_{2}(n,p)=\frac{\widehat{\overline{{s}}^{\prime}}({{r}}_{\text{c}}^{\prime },n,p){\text{I}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]+\frac{{\overline{{Q}}}_{1}^{\prime}(p)A(n){\text{I}}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]}{2\text{π}{{r}}_{\text{w}}^{\prime }(l-d){K}_{\text{r,0}}{\xi }_{0}(n,p)}}{{\text{I}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{K}}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]+{\text{K}}_{1}\left[{\xi }_{0}(n,p){{r}}_{\text{w}}^{\prime }\right]{I}_{0}\left[{\xi }_{0}(n,p){{r}}_{\text{c}}^{\prime }\right]}。$

(25c)

式中：I_j(·) (j = 0, 1)和K_j(·) (j = 0, 1)分别为第一类和第二类j阶虚宗量Bessel函数。对式（25）作有限Fourier余弦逆变换，则定解问题在Laplace空间的解可写作

$\bar s^{\prime}(r,z,p) = \frac{1}{B}\hat{\bar{s}^{\prime}}(r,0,p) + \frac{2}{B}\sum\limits_{n = 1}^\infty {\hat{\bar{s}^{\prime}}(r,n,p)} \cdot \cos ({w_n}z)\;\;。$

(26)

式中， ${\bar Q^{\prime}_1}(p)$ 的表达式仍然未知，可通过将式（25）和式（26）代入边界条件式（22）确定其表达式。需要注意的是，在计算 ${\bar Q^{\prime}_1}(p)$ 时还需首先确定受堵塞区外部边界的水头抬升 ${s_{\text{c}}}(z,t)$ ， ${s_{\text{c}}}(z,t)$ 又取决于通过过滤器的回灌流量Q(t)，而Q(t)的表达式仍然未知。为此，可以首先指定Q(t)一个初始值并计算该初始值条件下的 ${s_{\text{c}}}(z,t)$ ，进一步根据式（26）和 ${\bar Q^{\prime}_1}(p)$ 的表达式对Q(t)进行更新，不断循环迭代即可最终确定Q(t)的近似表达式。在迭代过程中，Q(t)可用下式近似：

$Q(t){\text{ = }}{Q_\infty }(t) + \left[ {{Q_0}(t) - {Q_\infty }(t)} \right] \cdot \sum {\left[ {{B_i}\exp ( - {\lambda _i}t)} \right]} 。$

(27)

式中，Q₀(t)和Q_∞(t)分别为井周含水层渗透系数为K_{r, 0}和K_{r, ∞}时的回灌流量，其表达式与Yang等^[3]推导的承压含水层中非完整定降深两区模型的抽水流量解相同，B_i和 ${\lambda _i}$ 为待定拟合参数，式（27）在i = 1且B_i = 1， ${\lambda _1}$ = ${\lambda _0}$ 时的计算结果可作为计算Q(t)的迭代初始值。

1.4 回灌初期与后期的近似解

在回灌伊始堵塞效应尚不明显，井周含水层与主含水层的导水特性无明显差异。此时，描述定水头回灌条件下井周含水层渗流场的定解问题与同等条件下主含水层渗流场的定解问题完全相同。将定解问题（2）~（9）中的K_r(t)和K_z(t)分别用K_{r, 0}和K_{z, 0}替换，则式（11）可以简化为 $t^{\prime} = t$ 。同时，令r_c → ∞，无穷远处的水头抬升可以忽略不计，即s_c(z, t) = 0。保持其他条件不变，得到新的定解问题。重复前述变量代换以及积分变换过程，得到新定解问题在Laplace空间的解为

$\begin{array}{l} \bar s(r^{\prime},z,p) = \frac{1}{B}\frac{{\bar q^{\prime}(p)A(0){{\text{K}}_0}\left[ {{\xi _0}(0,p)r^{\prime}} \right]}}{{{\xi _0}(0,p){K_{{\text{r,0}}}}{{\text{K}}_1}\left[ {{\xi _0}(0,p){{r}_{\text{w}}^{\prime}}} \right]}} + \\ \;\;\;\;\;\;\;\;\;\;\frac{2}{B}\sum\limits_{n = 1}^\infty {\frac{{\bar q^{\prime}(p){{\text{K}}_0}\left[ {{\xi _0}(n,p)r^{\prime}} \right]A(n)}}{{{\xi _0}(n,p){K_{{\text{r,0}}}}{{\text{K}}_1}\left[ {{\xi _0}(n,p){{r}_{\text{w}}^{\prime}}} \right]}}} \cos ({w_n}z) \;\;。 \end{array}$

(28a)

式中， $\bar q^{\prime}(p) = {{\bar Q^{\prime}(p)} \mathord{\left/ {\vphantom {{\bar Q^{\prime}(p)} {\left[ {2{\text{π }}{{r}_{\text{w}}^{\prime}}(l - d)} \right]}}} \right. } {\left[ {2{\text{π }}{{r}_{\text{w}}^{\prime}}(l - d)} \right]}}$ ，可根据边界条件式（22）确定，其表达式为

$\begin{aligned} \bar{q}(p)= & \frac{s_{\mathrm{w}}}{p}\left\{\frac{1}{B} \frac{\mathrm{K}_0\left[\xi_0(0, p) r_{\mathrm{w}}^{\prime}\right]}{\xi_0(0, p) K_{\mathrm{r}, 0} \mathrm{~K}_1\left[\xi_0(0, p) r_{\mathrm{w}}^{\prime}\right]}+\right. \\ & \left.\frac{2}{B(l-d)} \sum\limits_{n=1}^{\infty} \frac{\mathrm{K}_0\left[\xi_0(n, p) r_{\mathrm{w}}^{\prime}\right][A(n)]^2}{\xi_0(n, p) K_{\mathrm{r}, 0} \mathrm{K}_1\left[\xi_0(n, p) r_{\mathrm{w}}^{\prime}\right]}\right\}^{-1}\;\;。 \end{aligned}$

(28b)

式（28）即为非完整井定水头回灌初期水头抬升在Laplace空间的解，也是同等回灌条件下不考虑堵塞效应的承压含水层水头抬升在Laplace空间的解，且与Feng等^[1]推导的变流量条件下承压含水层中水头降深公式相同。特殊地，该解在 $r^{\prime}{\text{ = }}{r^{\prime}_{\text{c}}}$ 处的表达式即为式（9）右端项 ${s^{\prime}_{\text{c}}}(z,t^{\prime})$ 在Laplace空间的表达式。

当回灌试验进行足够长时间，井周堵塞影响下的含水层导水特性早已稳定。井周含水层水平与竖向渗透系数分别等于K_{r, ∞}和K_{z, ∞}。将式（2）~（9）中的K_r(t)和K_z(t)分别用K_{r, ∞}和K_{z, ∞}替换，得到新的定解问题。根据式（10）进行关于r的坐标变换并重复前述Laplace变换和Fourier变换，得到新定解问题在Laplace空间的解，其待定系数C₁(n, p)和C₂(n, p)表示如下：

${C}_{1}(n,p)=\frac{\widehat{\overline{{s}}^{\prime}}({{r}}_{\text{c}}^{\prime },n,p){\text{K}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]-\frac{\overline{{Q}}^{\prime}(p)A(n){\text{K}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]}{2\text{π}{{r}}_{\text{w}}^{\prime }(l-d){K}_{\text{r,}\infty }{\xi }_{\infty }(n,p)}}{{\text{I}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{K}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]+{\text{K}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{I}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]}\text{，}$

(29a)

${C}_{2}(n,p)=\frac{\widehat{\overline{{s}}^{\prime}}({{r}}_{\text{c}}^{\prime },n,p){\text{I}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]+\frac{\overline{{Q}}^{\prime}(p)A\text{(}n{\text{)I}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]}{2\text{π}{{r}}_{\text{w}}^{\prime }(l-d){K}_{\text{r,}\infty }{\xi }_{\infty }(n,p)}}{{\text{I}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{K}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]+{\text{K}}_{1}\left[{\xi }_{\infty }(n,p){{r}}_{\text{w}}^{\prime }\right]{\text{I}}_{0}\left[{\xi }_{\infty }(n,p){{r}}_{\text{c}}^{\prime }\right]}。$

(29b)

式（25a），（26），（29a），（29b）即为非完整井定水头回灌后期水头抬升在Laplace空间的解，其中 $\hat{\bar{s}^{\prime}}({r^{\prime}_{\text{c}}},n,p)$ 的表达式与式（28a），（28b）相同， ${\xi _\infty }{(n,p)^2} = w_n^2 +$ $Sp/(B{K_{{\text{z}},\infty }})$ 。此时， $\hat{\bar{s}^{\prime}}({r^{\prime}_{\text{c}}},n,p)$ 表达式中的 $\bar Q^{\prime}(p)$ 仍然未知，可联立式（25a）~（25c），（26），（28a），（28b），（29a），（29b）确定其表达式。

由于上述Laplace空间的解形式复杂，无法直接利用解析方法求得其在实时空间的表达式。本文在对比多种Laplace变换数值反演方法后，选用兼具计算效率和精度的Stehfest^[12]数值反演方法计算实时空间的水头抬升。需要指出的是，Stehfest算法中求和级数的项数N对反演结果具有显著影响，试算结果表明，当8≤N≤16时不同N值对应的反演结果差别不大，因此本文在进行Laplace变换的数值反演假设N = 16。

2. 结果验证与分析

目前由于缺少相关测试数据，本节首先通过数值计算方法对本文推导的半解析解进行对比验证，随后重点分析非完整井定水头回灌过程中堵塞诱发的渗透衰减对含水层渗流响应的影响规律。本节用于的验证和分析的非完整井定水头回灌井流模型如图 1所示。如无特殊声明，本节在进行模型计算时选取的试井参数和含水层参数分别如表 1，2所示。

表 1 井身结构参数

Table 1. Parameters of well structure 单位: m

试井类型	回灌压力	井径	滤管长度	滤管顶端高度	滤管底端高度
回灌井	200	0.1	10	15	5

下载: 导出CSV

| 显示表格

表 2 含水层水文地质参数

Table 2. Hydrological parameters of recharged aquifer

堵塞区厚度/m	初始渗透系数/(m·h^-1)		渐近渗透系数/(m·h^-1)		渗透衰减速率/h^-1	含水层厚度/m	储水系数
堵塞区厚度/m	水平	竖向	水平	竖向	渗透衰减速率/h^-1	含水层厚度/m	储水系数
1.0	0.06	0.006	0.03	0.003	0.05	20	0.004

下载: 导出CSV

| 显示表格

2.1 解的验证

当l = B和d = 0时，本文建立的定解问题退化为完整井定水头回灌井流模型。此时，含水层中的渗流为径向一维流动，回灌引起的水头抬升与深度无关。这时如果对Laplace域内堵塞区外部水头抬升 $\bar s^{\prime}({r^{\prime}_{\text{c}}},p)$ 作有限Fourier余弦变化，则可得到

$\hat{\bar{s}^{\prime}}\left(r_{\mathrm{c}}^{\prime}, n, p\right)=\left\{\begin{array}{ll} B \cdot \bar{s}^{\prime}\left(r_{\mathrm{c}}^{\prime}, p\right) & (n=0) \\ 0 & (n \neq 0) \end{array}\right. \;\;。$

(30)

此外，将l = B，d = 0代入A(n)化简得到

$A(n)=\left\{\begin{array}{ll} B & n=0 \\ 0 & n \neq 0 \end{array}\right. \;\;。$

(31)

将式（30），（31）代入式（25），（26）即得到完整井定水头回灌井流模型在Laplace空间的半解析解。此外，完整井定水头回灌井流模型还可通过MATLAB的偏微分方程求解函数pdepe进行数值求解。该函数可以求解如下形式的偏微分定解问题：

$c\left(x,t,u,\frac{\partial u}{\partial x}\right)\frac{\partial u}{\partial t}=\frac{1}{x}\frac{\partial }{\partial x}\left(xf\left(x,t,u,\frac{\partial u}{\partial x}\right)\right)+s\left(x,t,u,\frac{\partial u}{\partial x}\right)\text{，}$

(32a)

$u(x,0) = {u_0}(x) \;\;\text{，}$

(32b)

$p(x,t,u) + q(x,t)f(x,t,u,{{\partial u} \mathord{\left/ {\vphantom {{\partial u} {\partial x}}} \right. } {\partial x}}) = 0 \;\;。$

(32c)

其中，式（32a）~（32c）分别为定解问题的控制方程、初值条件和边值条件，x和t分别为空间变量和时间变量，c(x, t, u, ∂u/∂t)为控制方程非定常项的系数矩阵，f(x, t, u, ∂u/∂t)和s(x, t, u, ∂u/∂t)分别为控制方程的流项和源项。将完整井定水头回灌井流模型按照式（32a）~（34c）进行改写，并编写MATLAB计算程序，得到定解问题的数值解。图 2展示了不同位置的水头抬升的半解析计算值与数值计算值的对比结果。对比结果表明，利用Stehfest算法反演得到的实时空间的解在N = 16时与相应数值解符合良好，证明本文推导半解析解的准确性以及Stehfest算法在解决考虑变渗透系数的定水头非稳定流问题的适用性。

图 2 本文半解析解与数值解的比较

Figure 2. Comparison between semi-analytical and numerical solutions

下载: 全尺寸图片幻灯片

2.2 回灌堵塞程度的影响

图 3所示为不同堵塞程度，即不同K_{r, ∞}值，对井周堵塞区内不同位置（z = 10 m，r = 0.2，0.5，1.0 m）处水头抬升的影响，其中图 3（a）为水头抬升s时程曲线，图 3（b）为考虑渗透衰减（K_{r, ∞} < 0.06 m/h）与不考虑渗透衰减（K_{r, ∞} = 0.06 m/h）情形下的水头差Δs时程曲线。在回灌初期，堵塞效应尚不明显，任一给定位置处不同K_{r, ∞}值下的s随时间以相同速度逐渐增大。在该阶段不同位置处考虑堵塞效应与否的水头差Δs基本为零。随着回灌的进行，K_{r, ∞} = 0.06 m/h时的s继续随时间逐渐增大，而堵塞效应下K_{r, ∞} < 0.06 m/h时的s则呈现一定下降趋势，且K_{r, ∞}越小，r越大，下降趋势越显著，Δs越大。至回灌后期，堵塞引起的渗透衰减趋于稳定，此时不同K_{r, ∞}下的s均随时间缓慢增大，Δs在达到峰值后基本保持不变或轻微下降。

图 3 不同K_{r, ∞}值下堵塞区内水头抬升s和水头差

$\Delta s$ 变化规律

Figure 3. Development of s and

$\Delta s$ under different values of K_{r, ∞}

下载: 全尺寸图片幻灯片

图 4所示为不同K_{r, ∞}值对含水层水头抬升分布的影响，其中图 4（a）为不同时刻z = 10 m深度处的水头抬升s的径向分布曲线。从图 4中可以看出，任一给定时刻t的水头抬升随径向距离r逐渐减小，且K_{r, ∞}越小，s下降趋势越显著，s的大小及分布范围也越小。同时，水头抬升在沿径向距离下降的过程中，在堵塞区域与主含水层交界面处存在明显的转折，K_{r, ∞}越小，t越大，该转折越明显。对考虑渗透衰减与否两种情形下所得水头抬升分布曲线作差，计算结果如图 4（b）所示。可以发现，水头差分布沿径向距离先迅速增大，在井周含水层与主含水层交界面达到峰值，且K_{r, ∞}越小，t越大，水头差峰值越大，之后随着r的增大逐渐下降至零。

图 4 不同K_{r, ∞}值下含水层水头抬升s和水头差

$\Delta s$ 分布规律

Figure 4. Distribution of s and

$\Delta s$ under various values of K_{r, ∞}

下载: 全尺寸图片幻灯片

图 5所示为不同K_{r, ∞}对定水头条件下试验井回灌流量的影响，其中为回灌流量Q时程曲线，为考虑渗透衰减与否所得回灌流量差 $\Delta Q$ 时程曲线。从图 5中可以看出，堵塞诱发的渗透衰减效应显著降低了回灌流量。在回灌初期，渗透衰减效应尚不明显，不同K_{r, ∞}下回灌流量相差非常微小；随着渗透衰减效应的逐渐显现，K_{r, ∞}越小，回灌流量衰减越迅速，且考虑渗透衰减与否引起的流量差随着时间增加明显；当回灌时间足够长，渗透系数不再继续衰减，含水层中的渗流趋于拟稳态，此时不同K_{r, ∞}下的Q随着时间非常缓慢地下降，而 $\Delta Q$ 也在非常缓慢地下降。

图 5 不同K_{r, ∞}值下回灌流量Q和流量差

$\Delta Q$ 变化规律

Figure 5. Development of Q and

$\Delta Q$ under different values of K_{r, ∞}

下载: 全尺寸图片幻灯片

2.3 回灌堵塞速率的影响

所示为不同堵塞速率，即不同 $\lambda$ 值，对z = 10 m深度处井周堵塞区内含水层水头抬升的影响，其中图 6（a）为不同径向距离（r = 0.2，0.5，1.0 m）处水头抬升s时程曲线。从中可以看出，任一给定位置处的s在发展过程均存在显著的下降阶段，且 $\lambda$ 越大，r越大，下降阶段发生越早，下降程度越显著。在回灌初期和后期，不同的 $\lambda$ 值对含水层水头抬升几无影响，而在回灌中期，较大的 $\lambda$ 值可降低水头抬升。将 $\lambda$ ≠ 0.06/h时与 $\lambda$ = 0.06/h时的s作差，计算结果如所示。可以发现，相对于 $\lambda$ = 0.06/h时的水头抬升，越小的 $\lambda$ 引起的水头差Δs随时间增加越迅速、所能达到峰值越大、达到峰值所需时间越长。同时，Δs峰值大小和达到峰值所需的时间分别与r成正比和反比关系。此后，随着不同 $\lambda$ 下的渗透系数趋于相同的K_{r, ∞}，由于 $\lambda$ 不同引起的Δs随时间逐渐衰减至零。

图 6 不同

$\lambda$ 值下堵塞区内水头抬升s和水头差

$\Delta s$ 变化规律

Figure 6. Development of s and

$\Delta s$ under different values of

$\lambda$

下载: 全尺寸图片幻灯片

所示为井周含水层不同 $\lambda$ 值对不同时刻含水层水头抬升s分布的影响，其中为s的径向分布曲线，为因 $\lambda$ 不同而引起的水头抬升差Δs的径向分布曲线。从图中可以看出，任一给定时刻t的水头抬升随径向距离r逐渐减小。在回灌初期（t = 10 h）和后期（t = 1000 h），不同 $\lambda$ 下的s分布曲线接近重合在一起，此时的Δs分布几乎等于零。而在回灌中期（t = 50 h）， $\lambda$ 越小，任意r处的s越大，相应的s分布范围也越大。此时， $\lambda$ 不同引起的Δs随r先增大并在井周含水层与主含水层交界面达到最大值，且 $\lambda$ 越小，Δs随r增速越大、所能达到的峰值越大。此后，随着r继续增大，Δs开始缓慢下降并逐渐减少至零。

图 7 不同

$\lambda$ 值下含水层水头抬升s和水头差

$\Delta s$ 分布规律

Figure 7. Distribution of s and

$\Delta s$ under various values of

$\lambda$

下载: 全尺寸图片幻灯片

所示为不同 $\lambda$ 对定水头条件下试验井回灌流量的影响，其中为不同 $\lambda$ 情形下的回灌流量时程曲线。从图中可以看出， $\lambda$ 的大小对于试验初期和后期的回灌量几乎没有影响，但在试验中期，较大的 $\lambda$ 可一定程度上降低回灌流量。将 $\lambda$ ≠ 0.06/h时与 $\lambda$ = 0.06/h的回灌流量时程曲线作差，计算结果如所示。因 $\lambda$ 不同而引起的回灌流量差首先随着时间增大， $\lambda$ 越小，流量差的增加速度越快，流量差幅值越大，达到峰值所需的时间越短。此后，随着试验的进行，不同 $\lambda$ 计算得到的流量差逐渐减小为零。

图 8 不同

$\lambda$ 值下回灌流量Q和流量差

$\Delta Q$ 变化规律

Figure 8. Development of Q and

$\Delta Q$ under different values of

$\lambda$

下载: 全尺寸图片幻灯片

3. 分析与讨论

本文的参数分析表明，堵塞程度同时改变了含水层水头抬升速度以及含水层的最终渗流状态，而堵塞速率则主要改变了含水层的渗流响应速度，并不影响含水层最初和最终的渗流状态。当井周含水层发生堵塞时，含水层渗透性下降增加了回灌水从试井流入含水层的难度，导致回灌流量的下降，进而引起整个含水层水头抬升不同程度的降低。同时，由于井筒定水头边界的约束作用，这种堵塞引起的水头抬升下降在井筒附近非常微小。随着远离井筒，水头抬升下降逐渐增大并在堵塞区与主含水层交界面达到最大值。此后，随着径向距离和过水断面的增大，井周堵塞导致的水头差逐渐减小并最终降为零。

上文的参数分析同时还表明，当考虑堵塞引起的渗透系数衰减时，含水层中水头抬升呈现先快速上升，继而明显下降，随后又继续缓慢上升的发展过程。这一现象与Wen等^[13]报道的抽水流量呈指数形式衰减时引起的水头变化过程非常类似。结合Wen等^[13]的研究，可以推断这种水头抬升过程中的下降趋势是由于渗透系数衰减导致恒定压力下回灌流量下降引起的。当堵塞越严重、回灌流量衰减越显著时，这种水头抬升的下降趋势越明显。由于该特征与渗透性恒定含水层中进行定水头回灌引起的承压水头持续抬升存在显著差异，这一特征可以作为判别井周含水层是否发生堵塞以及堵塞程度的判断标准。

除了水头抬升过程中出现的下降现象，井周含水层渗透性衰减还导致了水头抬升分布在井周堵塞区与主含水层交界面存在一定的折点，且该折点随着堵塞的发展愈加明显。分布曲线折点的存在主要是由于井周含水层与主含水层存在渗透差异造成的，这与均质含水层中的光滑分布曲线形成明显区别。因此，这一特征可以用于判别堵塞范围r_c的大小。需要指出的是，由于本文建立的井流模型假设r_c为定值，上述r_c的确定方法仅适用于试井周围堵塞分布在水平和竖直方向上相对均匀的情形，而对于实际工程中可能存在的堵塞分布不均匀的问题，相关模型的建立以及堵塞范围的判定还有待进一步研究。

4. 结论

本文以井周含水层为研究对象，建立了考虑井周堵塞效应的非完整井定水头回灌井流模型。通过变量代换、Laplace变换和有限余弦Fourier变换，推导了井流模型在Laplace-Fourier空间的解，利用有限余弦Fourier逆变换及基于Stehfest法的Laplace数值逆变换得到了井流模型在实时空间的解，并通过与数值解对比证明了本文求解方法的正确性。在此基础上，本文分析了非完整井定水头回灌条件下堵塞效应对渗流场的影响规律，得到以下3点结论。

（1）较小的K_{r, ∞}会减小回灌过程中水头抬升和回灌流量，并降低含水层的拟稳态渗流响应，而较大的 $\lambda$ 仅降低了回灌过程中水头抬升和回灌流量，并不影响含水层初始和达到拟稳态时的渗流状态。

（2）受定水头回灌边界的影响，因K_{r, ∞}和λ不同而引起的水头差沿径向距离先增大后减小，在井周堵塞区与主含水层交界面达到最大值，较小的K_{r, ∞}或较大的λ均会引起水头差分布与时程曲线的增大。

（3）回灌堵塞使得水头抬升分布在井周堵塞区与主含水层交界面上存在明显的折点，造成水头在抬升过程中存在明显的下降阶段，上述特征在堵塞越严重时越显著，可为回灌堵塞的评估与预测提供参考。

图 1 水布垭Ⅰ期面板裂缝分布

Figure 1. Cracks in Stage Ⅰ slab of Shuibuya CFRD

下载: 全尺寸图片幻灯片

图 2 HJD大坝左岸水平向裂缝分布

Figure 2. Horizontal cracks in left abutment slab of HJD CFRD

下载: 全尺寸图片幻灯片

图 3 JL面板坝面板裂缝分布

Figure 3. Cracks in slab of JL CFRD

下载: 全尺寸图片幻灯片

图 4 多点约束法三维示例

Figure 4. 3D example of multi-point constraint refinement

下载: 全尺寸图片幻灯片

图 5 钢筋埋置单元

Figure 5. Steel bar embedding element

下载: 全尺寸图片幻灯片

图 6 混凝土单轴应力应变曲线

Figure 6. Uniaxial stress-strain curve of concrete

下载: 全尺寸图片幻灯片

图 7 大石峡面板坝防渗体布置

Figure 7. Arrangement of impervious body of Dashixia CFRD

下载: 全尺寸图片幻灯片

图 8 大石峡大坝材料分区

Figure 8. Material zoning of Dashixia CFRD

下载: 全尺寸图片幻灯片

图 9 基于非协调网格的大石峡面板坝计算模型

Figure 9. Refined FEM model for Dashixia CFRD

下载: 全尺寸图片幻灯片

图 10 面板内钢筋分布

Figure 10. Reinforcements in slab

下载: 全尺寸图片幻灯片

图 11 蓄水运行期面板变形分布

Figure 11. Deformations of slab during operation period

下载: 全尺寸图片幻灯片

图 12 蓄水运行期面板应力分布

Figure 12. Slab stresses during operation period

下载: 全尺寸图片幻灯片

图 13 不同计算模型1640 m高程面板轴向应力比较

Figure 13. Comparison of axial stresses at elevation 1640 m by different models

下载: 全尺寸图片幻灯片

图 14 蓄水运行期面板受拉损伤参数d_t分布

Figure 14. Tensile damage factor d_t of slab during operation period

下载: 全尺寸图片幻灯片

图 15 钢筋拉应力分布

Figure 15. Tensile stresses of reinforcements

下载: 全尺寸图片幻灯片

图 16 面板裂缝分布

Figure 16. Cracks in slabs of Dashixia CFRD

下载: 全尺寸图片幻灯片

图 17 面板配筋加强区示意

Figure 17. Special reinforcement area of Dashixia CFRD

下载: 全尺寸图片幻灯片

表 1 筑坝料“南水”模型参数

Table 1 Parameters of "NHRI" model

坝体分区	${\rho _{\text{d}}}$ /(g·cm^-3)	${\varphi _0}$ /(°)	$\Delta \varphi$ /(°)	$k$	$n$	${R_{\text{f}}}$	c_d/%	n_d	R_d
垫层区（2A）	2.31	49.3	5.9	920	0.29	0.67	0.49	0.45	0.65
主砂砾石区（3BA、3BB）	2.27	50.1	6.3	1294	0.32	0.74	0.30	0.67	0.72
主堆石区（3BC、3C2）	2.22	53.2	9.0	1104	0.22	0.65	0.40	0.72	0.64
下游堆石区（3C1）	2.22	52.6	8.7	1044	0.21	0.66	0.46	0.70	0.66
增模区	2.28	69.1	17.9	0.14	3928.7	0.35	0.08	0.44	0.11

下载: 导出CSV

表 2 筑坝料流变模型参数

Table 2 Parameters of creep model

坝体分区	α	b /%	c/%	d /%	m₁	m₂	m₃
垫层区（2A）	0.006	0.130	0.028	0.277	0.355	0.600	0.700
主砂砾区（3BA、3BB）	0.0062	0.114	0.025	0.241	0.351	0.578	0.696
主堆石料区（3BC、3C2）	0.0055	0.138	0.030	0.363	0.358	0.629	0.783
下游堆石区（3C1）	0.0054	0.140	0.032	0.365	0.358	0.629	0.783

下载: 导出CSV

表 3 混凝土损伤模型参数

Table 3 Parameters for damage model

E_c/10⁴ MPa	f_cr/MPa	ε_cr	α_c	f_tr/MPa	ε_tr	α_t
3.0	20.1	1.470×10^-3	0.74	2.01	9.5×10^-5	1.25

下载: 导出CSV

表 4 不同配筋方案面板裂缝宽度

Table 4 Cracks width under different reinforcements scheme

	配筋加强区Ⅰ		配筋加强区Ⅱ
	水平向配筋率/%	纵向裂缝最大宽度/mm	顺坡向配筋率/%	水平向最大裂缝宽/mm
原方案	0.587	0.32	0.587	0.24
加强方案	0.783	0.20	0.691	0.18

下载: 导出CSV

参考文献(22)

[1]	罗福海, 张保军, 夏界平. 水布垭面板堆石坝施工期裂缝成因及处理措施[J]. 水利水电快报, 2010, 31(12): 5-8. https://www.cnki.com.cn/Article/CJFDTOTAL-SLSK201012002.htm LUO Fuhai, ZHANG Baojun, XIA Jieping. Causes and treatment measures of cracks in Shuibuya concrete face rockfill dam during construction[J]. Express Water Resources & Hydropower Information, 2010, 31(12): 5-8. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLSK201012002.htm
[2]	混凝土面板堆石坝设计规范: N/T 10871—2021[S]. 北京: 中国水利水电出版社, 2022. Code for Design of Concrete Face Rockfill Dams: N/T 10871—2021[S]. Beijing: China Water & Power Press, 2022. (in Chinese)
[3]	钟红, 林皋, 胡志强. 有限元计算中疏密网格过渡方法研究[J]. 计算力学学报, 2007, 24(6): 887-891, 898. https://www.cnki.com.cn/Article/CJFDTOTAL-JSJG200706030.htm ZHONG Hong, LIN Gao, HU Zhiqiang. Two methods for transition between coarse and fine finite elements[J]. Chinese Journal of Computational Mechanics, 2007, 24(6): 887-891, 898. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JSJG200706030.htm
[4]	周墨臻, 张丙印, 王伟. 高面板堆石坝软缝接触计算模型及其数值实现[J]. 岩石力学与工程学报, 2016, 35(增刊1): 2803-2810. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2016S1024.htm ZHOU Mozhen, ZHANG Bingyin, WANG Wei. Algorithm and simulation methods for the soft longitudinal joint of the concrete faces in high concrete-faced rockfill dams[J]. Chinese Journal of Rock Mechanics and Engineering, 2016, 35(S1): 2803-2810. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX2016S1024.htm
[5]	周墨臻, 张丙印, 张宗亮, 等. 超高面板堆石坝面板挤压破坏机理及数值模拟方法研究[J]. 岩土工程学报, 2015, 37(8): 1426-1432. doi: 10.11779/CJGE201508010 ZHOU Mozhen, ZHANG Bingyin, ZHANG Zongliang, 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) doi: 10.11779/CJGE201508010
[6]	孔宪京, 陈楷, 邹德高, 等. 一种高效的FE-PSBFE耦合方法及在岩土工程弹塑性分析中的应用[J]. 工程力学, 2018, 35(6): 6-14. https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201806004.htm KONG Xianjing, CHEN Kai, ZOU Degao, et al. An efficient fe-psbfe coupled method and its application to the elasto-plastic analysis of geotechnical engineering structures[J]. Engineering Mechanics, 2018, 35(6): 6-14. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCLX201806004.htm
[7]	邹德高, 刘锁, 陈楷, 等. 基于四叉树网格和多边形比例边界有限元方法的岩土工程非线性静动力分析[J]. 岩土力学, 2017, 38(增刊2): 33-40. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2017S2005.htm ZOU Degao, LIU Suo, CHEN Kai, et al. Nonlinear static and dynamic analysis for geotechnical engineering based on quadtree mesh and polygon scaled boundary finite element method[J]. Rock and Soil Mechanics, 2017, 38(S2): 33-40. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2017S2005.htm
[8]	邹德高, 陈楷, 刘锁, 等. 非线性比例边界有限元在面板坝分析中的应用[J]. 土木与环境工程学报, 2019, 41(3): 11-18. https://www.cnki.com.cn/Article/CJFDTOTAL-JIAN201903002.htm ZOU Degao, CHEN Kai, LIU Suo, et al. Application of nonlinear scaled boundary polygon element method in analysis of concrete face rockfill dam[J]. Journal of Civil and Environmental Engineering, 2019, 41(3): 11-18. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JIAN201903002.htm
[9]	魏匡民, 陈生水, 李国英, 等. 位移多点约束法在面板堆石坝精细模拟中的应用研究[J]. 岩土工程学报, 2020, 42(4): 616-623. doi: 10.11779/CJGE202004003 WEI Kuangmin, CHEN Shengshui, LI Guoying, et al. Application of displacement multi-point constraint refinement method in simulation of concrete-faced rockfill dams[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(4): 616-623. (in Chinese) doi: 10.11779/CJGE202004003
[10]	王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2003. WANG Xucheng. Finite Element Method[M]. Beijing: Tsinghua University Press, 2003. (in Chinese)
[11]	陈国荣. 有限单元法原理及应用[M]. 北京: 科学出版社, 2009. CHEN Guorong. Principle and Application of Finite Element Method[M]. Beijing: Science Press, 2009. (in Chinese)
[12]	孔宪京, 屈永倩, 邹德高, 等. 强震作用下面板堆石坝跨尺度面板开裂演化分析[J]. 岩土工程学报, 2020, 42(6): 989-996. doi: 10.11779/CJGE202006001 KONG Xianjing, QU Yongqian, ZOU Degao, et al. Cross-scale crack evolution analysis for face slab in concrete faced rockfill dams under strong earthquake[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(6): 989-996. (in Chinese) doi: 10.11779/CJGE202006001
[13]	巫昌海, 汪基伟. 混凝土三维钢筋埋置组合式有限单元模型及其网格自动生成[J]. 计算机辅助设计与图形学学报, 2000, 12(10): 761-764. https://www.cnki.com.cn/Article/CJFDTOTAL-JSJF200010011.htm WU Changhai, WANG Jiwei. Three-dimensional embedded finite element model of reinforced concrete and its automatic mesh generation[J]. Jouraal of Computer Aided Design & Computer Graphics, 2000, 12(10): 761-764. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JSJF200010011.htm
[14]	XU B, ZOU D G, KONG X J, et al. Dynamic Damage evaluation on the slabs of the concrete faced rockfill dam with the plastic-damage model[J]. Computers and Geotechnics, 2015, 65: 258-265. http://www.onacademic.com/detail/journal_1000037433272510_9a3e.html
[15]	CEN W J, WEN L S, ZHANG Z Q, et al. Numerical simulation of seismic damage and cracking of concrete slabs of high concrete face rockfill dams[J]. Water Science and Engineering, 2016, 9(3): 205-211. http://qikan.cqvip.com/Qikan/Article/Detail?id=7000080633
[16]	混凝土结构设计规范: GB 50010—2010[S]. 北京: 中国建筑工业出版社, 2011. Code for Design of Concrete Structures: GB 50010—2010[S]. Beijing: China Architecture & Building Press, 2011. (in Chinese)
[17]	朱百里, 沈珠江. 计算土力学[M]. 上海: 上海科学技术出版社, 1990. ZHU Baili, SHEN Zhujiang. Computational Soil Mechanics[M]. Shanghai: Shanghai Science and Technology Press, 1990. (in Chinese)
[18]	李国英, 米占宽, 傅华, 等. 混凝土面板堆石坝堆石料流变特性试验研究[J]. 岩土力学, 2004, 25(11): 1712-1716. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200411008.htm LI Guoying, MI Zhankuan, 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
[19]	顾淦臣, 沈长松, 岑威钧. 土石坝地震工程学[M]. 北京: 中国水利水电出版社, 2009. GU Ganchen, SHEN Changsong, CEN Weijun. Earthquake Engineering for Earthrock Dams[M]. Beijing: China Water & Power Press, 2009. (in Chinese)
[20]	高瑞平. 混凝土构件裂缝宽度计算方法研究[D]. 南京: 东南大学, 2008. GAO Ruiping. Study on Calculation Method of Crack Width of Concrete Members[D]. Nanjing: Southeast University, 2008. (in Chinese)
[21]	丁旭柳, 伍鹤皋, 朱忠华. 钢筋混凝土压力隧洞裂缝宽度计算方法的探讨[J]. 武汉大学学报(工学版), 2001, 34(2): 24-27. https://www.cnki.com.cn/Article/CJFDTOTAL-WSDD200102006.htm DING Xuliu, WU Hegao, ZHU Zhonghua. Study on computational method of crack width of reinforced concrete pressure tunnels[J]. Engineering Journal of Wuhan University, 2001, 34(2): 24-27. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-WSDD200102006.htm
[22]	兰伟钦, 肖明. 考虑钢筋滑移效应的高压隧洞衬砌配筋计算方法[J]. 水电与新能源, 2021, 35(4): 66-71. https://www.cnki.com.cn/Article/CJFDTOTAL-HBFD202104016.htm LAN Weiqin, XIAO Ming. On the reinforcement calculation of high-pressure tunnel lining considering the slip effect of reinforcement[J]. Hydropower and New Energy, 2021, 35(4): 66-71. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-HBFD202104016.htm