抽降水引发的弱透水层非线性固结解析解

2021-06-10李传勋江留慧

工程科学与技术 2021年3期

李传勋，江留慧

(江苏大学土木工程与力学学院，江苏镇江 212013)

随着经济发展及城市化进程的加快，人类对地下水资源的需求日趋增长，其后果是由大规模抽取地下水引发的地面沉降问题日益严重[1–3]。大范围地面沉降给人类生活和生产环境、生命财产等造成危害，并严重制约社会和经济的可持续发展。因此，深入研究由抽取地下水引发的地面沉降问题具有重要的理论和实际意义。

计算机技术的发展极大地促进了Biot固结理论在抽降水引发地面沉降研究中的应用。Lewis等[4]最早建立基于Biot耦合固结理论的有限元模型,并分析了威尼斯地面沉降问题。Bear等[5]也基于Biot固结理论建立含水层抽水引起的地面沉降计算模型，分析了降水引起的区域地面沉降问题。Teatini等[6]利用基于Biot固结理论建立的3维有限元模型，分析了Emilia–Romagna地区的地面沉降问题。陈杰等[7]运用2维Biot固结平面有限元程序计算因超采地下水而引起的地面沉降问题。骆祖江等[8]以Biot固结理论为基础，考虑土体的非线性压缩及渗透特性，建立了浅层地下水开采与地面沉降的3维全耦合模型。由此可见，地下水位变化引起的地面沉降采用Biot 3维固结理论建立的完全耦合模型计算已相对完善。但也需注意到Biot 3维固结理论必须借助于相对复杂的数值计算，模型较复杂、参数较多导致数值计算不仅耗时，而且计算结果往往因人而异，很难在实际工程中得到推广应用。

事实上，抽水引起的场地地面沉降范围与发生主要沉降变形的土层厚度相比往往大几倍、几十倍之多，此时采用1维固结理论计算弱透水层沉降是完全适用的。鉴于此，不少学者考虑土体压缩性和渗透性参数与孔隙比的关系，建立了3维地下水渗流和1维固结理论组合的部分耦合模型，对降水引起的地面沉降进行数值模拟[9–11]，但仍不可避免数值计算的局限性。

解析解能简捷地预测地面沉降的时空发展规律，便于工程应用。故学者们相继开展了由抽降水引发的地面沉降相关解析解的研究。如：基于太沙基1维固结理论，骆冠勇等[12]推导了承压水大面积减压引起的上覆弱透水层1维线性固结解析解。谢康和等[13]获得了越流系统中承压层水位大面积瞬时下降引发的弱透水层的1维固结解析解。谢海澜等[14]考虑土体固结过程中的非达西渗流特性，给出了承压层水位瞬时下降的弱透水层固结半解析解。张玮鹏等[15]考虑土中渗流存在的起始比降，推导了承压层水位瞬时下降引发的弱透水层固结解析解。夏长青[16]考虑土体成层性，推导了潜水层和承压层水位随时间线性下降情况下成层软土固结解析解。晏凌晗[17]考虑土体流变特性，推导了潜水层水位瞬时下降引发的弱透水层1维流变固结解析解。上述研究虽然考虑了软黏土的成层性、流变性、非达西渗流特性等，但均假定土体在固结过程中其渗透系数和压缩系数保持不变，而实际土体在固结过程中土体渗透系数和压缩系数随有效应力增加而不断变化。

Li[18]基于非线性应力应变关系推导了越流系统中，水位波动变化时，软土层1维非线性固结解析解，但其未能考虑土体渗透性的变化。张云等[19]考虑了土体的非线性变形及渗透特性，建立了以有效应力为变量的1维非线性地面沉降模型，用半解析法进行了求解。黄大中[20]在压缩系数与渗透系数同比例减小（即固结系数保持不变）的假定下得到了潜水层和承压层水位随时间大面积线性下降时的1维非线性固结解析解。但实际上土体在固结过程中压缩系数与渗透系数的变化往往并不成比例。目前，对由抽降水引发的、能反映土体非线性固结特性的弱透水层非线性固结解析解的研究仍较缺乏。

综上所述，十分有必要深入开展抽降水条件下，考虑土体非线性压缩和渗透特性的固结理论研究，尤其是更适合工程技术人员使用的、形式简单的非线性固结解析解。本文引入土体经典的e–lg σ′和e–lgkv非线性关系描述土体的压缩和渗透特性，建立潜水层大面积抽水引发的下卧弱透水层的非线性固结模型，并获得其近似解析解。该解答不仅进一步完善了非线性固结理论体系，而且其简单易用，便于在实际工程中推广应用。

1 问题描述及固结控制方程建立

1.1 问题描述

图1为潜水层水位下降引发的弱透水层1维固结示意图。如图1所示：z为以弱透水层顶面为坐标原点、向下为正的竖向坐标系，弱透水层厚度为H，底面不透水；假定弱透水层土体初始有效应力沿深度保持不变，其值等于自重应力沿深度的加权平均值，将其记为σ′0；潜水层发生大面积抽水（降水），h0为t=0时刻初始变化水位，最终水位下降值为hc；潜水层水位下降引起弱透水层顶面处总应力和孔压均减小，但孔隙水压力减小更明显，故引发弱透水层中土体有效应力进一步增大，进而引发弱透水土层的固结沉降变形。

图1 潜水层水位下降引发的弱透水层1维固结示意图Fig. 1 Schematic diagram of one-dimensional consolidation of aquitard caused by the dropping of groundwater table in phreatic aquifer

图2为潜水层中水位变化h(t)与时间t的关系曲线。特殊的水位瞬时下降和单级等速下降模式如虚线所示。

图2 水位下降量与时间的关系Fig. 2 Relationship between the dropping of groundwater table and time

水位瞬时下降其表达式为：

水位单级等速下降时表达式为：

式中，tc为水位下降历时。

1.2 固结相关土工参数的推导

引入经典的e–lg σ′和e–lgkv关系描述土体固结过程中的非线性压缩与渗透特性：

式中，Cc、Ck分别为压缩指数和渗透指数，e为孔隙比，e0为初始孔隙比，σ′为土层竖向有效应力，kv为土体的渗透参数，和kv0分别为e–lg σ′和e–lgkv曲线上相对应于e0的土体初始有效应力和初始竖向渗透系数。

对式（3）求导，可得土体体积压缩系数mv为：

式（4）减式（3），可得土体渗透系数kv为：

进而根据土体固结系数Cv定义得到：

1.3 固结控制方程及求解条件

饱和土体1维非线性固结普遍连续方程为：

式中，u为超静孔隙水压力。

将式（3）对t求偏导可得：

根据有效应力原理，土体有效应力表达为：

式中，q(t)为外界因素改变引起的总应力增量。

将式（7）、（9）和（10）代入式（8）可得1维非线性固结控制方程为：

潜水层水位下降，计算弱透水层顶面处总应力时，水位下降区域内土体重度由饱和重度降为自然重度，导致作用在弱透水层顶面处总应力降低，t时刻总应力增量q为：

式中，γ为降水后潜水层中砂土的天然重度，γsat为潜水层中土的饱和重度。

弱透水层顶面处水头变化与潜水层水位变化一致，土层的边界条件表达为：

其初始条件表达为：

黄大中[20]在Cc/Ck=1（固结系数在固结过程中保持不变）下，获得该模型解析解。目前，关于Cc/Ck≠1时，潜水层水位下降引发的弱透水层1维非线性固结解析解研究甚少，本文将对此问题展开探究。

2 固结模型的近似解析解

2.1 变量代换后的控制方程及求解条件

式中，p0=γwh0。

2.2 超静孔隙水压力解答

根据控制方程及边界条件形式，假定满足控制方程和边界条件的解的形式为：

式中：M=(2m–1)π/2，m=1，2，3，…；βm、Bm、Cm为待定参数；Tm(t)为时间t的函数，由R(t)确定。

式中，τ为积分变量。

进而可得超静孔隙水压力的普遍解答为：

在式（23）基础上，可求得土层按孔压定义的平均固结度Up为：

t时刻土层发生的沉降变形值St为：

土层的最终沉降变形值S∞为：

进而得到按变形定义的土层平均固结度Us为：

3 特殊降水模式下的解答

应用以上变量，边界条件和初始条件分别转化为：

3.1 水位瞬时下降模式

图2中虚水平线代表t>0时，潜水层水位发生瞬时下降，水位下降深度为hc。可得该降水模式下总应力

增量和求解条件分别为：

将式（29）代入待定参数Bm和式（22），可解得：

将Bm、Cm、βm代入式（23）得到弱透水层中超静孔隙水压力为：

式中，Tv=Cv0t/H2。

3.2 水位单级等速下降模式

如图2中虚斜线所示，当t≤tc时，潜水层水位发生等速下降，t>tc时水位下降值恒定为hc。可得该降水方式下总应力增量和求解条件分别为：

同理，可求得超静孔压表达式为：

式中，Tvc=Cv0tc/H2。

3.3 Cc/Ck→1条件下超静孔压解的退化

在Cc/Ck→1的条件下，水位单级等速下降时的超静孔压解可退化为文献[20]中相同条件下Cc/Ck=1时解的形式。超静孔压解u的退化过程为：

4 工程实例对比

为验证本文解答的正确性，将其与某工程降水引发的地面沉降观测结果、分层总和法计算结果及数值模拟结果[22]进行对比。

简化文献[22]中第五章表5–2和图5–22所示工程实例土层信息，将杂填土、粉土、黏土、粉砂和细砂土层视为土质均一的潜水层，天然重度取各土层天然重度的厚度加权平均值为19.637 kN/m3，这种简化对于弱透水层的沉降计算无影响。初始水位为5.5 m；隔水层粉质黏土视为弱透水层，弱透水层底面不透水；土层性质详见表1，降水引起的有效应力增量如表2所示。在距离基坑5.2 m处，水位下降6.4 m，忽略潜水层水位下降区域内饱和土体重度变化，根据表2和式（26）可分别计算潜水层和弱透水层沉降量。

将计算结果与文献[22]对比，如表3所示。可以看出，本文解与分层总和法计算结果相比较更趋近于实测沉降值，可能原因主要体现两个方面：其一，本文认为水位下降区域最终有效应力增加值呈线性分布，其更能体现土体中的真实受力特性，故本文相较于分层总和法中有效应力保持恒值，更能体现土体的真实状况；其二，本文解弱透水层沉降计算考虑了土体压缩性与初始有效应力间的非线性关系。本文计算的沉降充分考虑土层中的实际有效应力，其结果更能反映土体的真实性状。从计算的总沉降量与实测结果对比说明了本文的弱透水层沉降计算解答是合理有效的。

表1 土性参数Tab. 1 Parameters of soil characteristic

表2 有效应力增量计算Tab. 2 Calculation of effective stress increment

表3 不同计算方法下沉降量对比Tab. 3 Comparison of settlement results with different calculation methods

5 算例及固结性状分析

某越流系统，砂土层饱和重度γsat=20 kN/m3，自然重度γ=18 kN/m3，水的重度按γw=10 kN/m3，弱透水层土层厚度H=10 m，顶面透水，底面不透水，弱透水层自重应力的平均值50 kPa。依据郑刚等[23]反演得到的弱透水层土性参数计算表，取弱透水层压缩指数Cc=0.02，土体初始孔隙比e0=0.7，与初始孔隙比e0相对应的渗透系数kv0=4×10–9m/s。Tv<0.05时潜水层水位等速下降，而后水位下降值恒定为hc=5 m。

基于上述算例，利用本文解分别计算绘制在不同Cc/Ck值、水位下降速率Tvc、水位下降终值hc及砂土层自然重度γ下，弱透水层固结沉降曲线，分析潜水层水位单级等速下降时弱透水层的1维非线性固结性状。

5.1 Cc/Ck对固结性状的影响

图3为不同Cc/Ck下，弱透水层平均固结度与时间因子的关系曲线。图4和5分别为改变Cc和Ck值时弱透水层沉降量与时间因子的关系曲线。如图3～5所示：Cc/Ck越大，相同时间因子下弱透水层的平均固结度越小，即固结速率随Cc/Ck增大而减慢；当Ck值保持不变时，Cc值越大，Cc/Ck越大，虽然此时固结速率会减慢，但相同时间因子下弱透水层的沉降量及最终沉降量仍随Cc/Ck增大而变大；当Cc值保持不变时，Ck值越小，Cc/Ck越大，固结完成前，相同时间因子下土层沉降量越小，但固结完成时的最终沉降量相同。该固结性状可由固结系数表达式（7）和沉降表达式（25）得到解释：Cc/Ck越大，固结系数越小，固结速率减慢，最终沉降量与Cc值成正比，与Ck值无关；Cc值越大，Cc/Ck越大，固结速率减慢，但相同时间因子下沉降量越大；当其他参数一定时，增大Ck，固结速率及沉降速率加快，但不改变沉降量终值。

图3 Cc/Ck对平均固结度Up曲线的影响Fig. 3 Influence of Cc/Ck on the curves of average consolidation degree Up

图4 Cc值对沉降曲线的影响Fig. 4 Influence of Cc on settlement curves

图5 Ck值对沉降曲线的影响Fig. 5 Influence of Ck on settlement curves

5.2 水位下降速率对固结性状的影响

图6为水位下降速率与沉降量的关系曲线。由图6可知：随着Tvc增大（水位下降速率减慢），同一时间因子下沉降量随之减小，但最终沉降量相同；在水位下降历时发生改变而其他参数不变的条件下，水位瞬时下降时弱透水层沉降速率最快。该固结性状可由沉降量表达式（25）解释，水位下降越快，在相同时间因子下，弱透水层有效应力越大，相应的沉降量越大。

图6 水位下降速率对沉降曲线的影响Fig. 6 Influence of the dropping rate of groundwater table on settlement curves

5.3 不同Cc/Ck下水位下降终值对固结性状的影响

图7、8分别为不同水位下降终值时，σ′/σ′0沿深度分布和沉降量随时间变化曲线。由图7～8可得：相同Cc/Ck和同一深度处的有效应力σ′随hc增大而增大；hc越小，Cc/Ck值对有效应力σ′的影响越小。相同Cc/Ck和时间因子下，hc越大，沉降量越大；相同hc和时间因子下，Ck越大，Cc/Ck越小，沉降量越小，但最终沉降量相同；Cc/Ck对沉降曲线影响随hc减小而减弱，且hc=1不同，Cc/Ck下沉降曲线几乎重合，甚至可忽略。由式（7）和（25）可知，Cc/Ck越小，固结系数越大，固结速率跟沉降速率越快；hc越大，相同时间因子下，有效应力越大，相应的沉降量越大；hc越小，σ′/σ′0越趋近于1，Cc/Ck对固结速率的影响越小，同时Ck不影响沉降量终值，故而出现上述性状。

图7 Cc/Ck下水位下降终值对σ′/σ′0沿深度分布的影响Fig. 7 Influence of the final dropping values of groundwater table on distribution of σ ′/σ′0 along depth with different Cc/Ck

图8 Cc/Ck下水位下降终值对沉降曲线的影响Fig. 8 Influence of the final dropping values of groundwater table on settlement curves with different Cc/Ck

5.4 砂土层自然重度γ对固结速率的影响

图9、10分别为不同砂土自然重度γ下，弱透水层超静孔压随时间和深度变化曲线。

图9 砂土层自然重度γ对超静孔压消散的影响Fig. 9 Influence of specific density of sand layer on dissipation of excess pore pressure

图 10 砂土层自然重度对超静孔压沿深度分布的影响Fig. 10 Influence of specific density of sand layer on distribution of excess pore pressure with depth

从图9和10可发现：相同深度处的超静孔隙水压力随γ增大而增大，说明孔压消散速率随γ增大而减慢；相同条件下，相较于Cc/Ck>1，Cc/Ck<1时，对应的超静孔压值更小，即Cc/Ck越小，孔压消散速率越快。

图11、12分别为不同砂土层自然重度γ下，Cv/Cv0沿深度分布曲线和沉降量随时间因子变化曲线。由式（7）和式（25）可知：Cc/Ck<1时，砂土层自然重度γ越大，相同时间因子下有效应力σ′越大，固结系数越大，固结速率越快，相应的沉降量越大；Cc/Ck>1时，砂土层自然重度γ越大，相同时间因子下有效应力σ′越大，固结系数越小，固结速率越慢，但相应的沉降量越大；γ越大，Nq越大，最终沉降量越大。