沉井施工降水时地下水动态分析与对江堤影响的研究毕业论文

沉井施工降水时地下水动态分析与对江堤影响的研究1导言1.1地下水模拟的重要性地下水动态分析实际上是地下水系统的数值模拟。地下水是宝贵的自然资源，是维持人类社会可持续发展的经济和战略资源，是综合国力的有机组成部分[1]。当今世界面临的“人口、资源、环境”三大问题都直接或间接地与地下水有关。地下水储水量大，水量稳定，不易被污染，水质好。现在江河湖泊污染严重，地下水在我国城市生活和工农业用水中占有重要地位。中国约有三分之二的城市将地下水作为重要的供水水源。随着人口增长、工农业生产和城市建设的快速发展，地下水资源的开发利用规模日益增大，由此产生了水质恶化、形成区域性降落漏斗、地面沉降等一系列地质灾害和环境问题，给人类带来了不同程度的危害。因此，有必要查明地下水资源的形成条件，做出正确的总量评价，对研究区地下水系统进行模拟分析，并根据其时空分布特征制定合理的开发利用方案[2]。地下水系统[3]是指在一定水文地质条件下，在一定围隔内形成的，水力联系紧密，与相邻含水系统相对隔绝的地下含水系统。从广义上讲，地下水系统是指包括地下含水层系统及其相关的社会、经济和环境因素的整体。实际水文地质条件复杂:地形起伏、含水层厚度不均、非均质、各向异性、结构复杂、多层不连续、源汇项时空分布不均、水文地质参数时空变异性明显等。数值模型对这些特征有相对合理的表示，是对真实地下水系统的模拟和仿真。地下水资源评价是“对一定系统中地下水资源的时空分布特征和开发利用条件作出科学全面的分析和估计”[4][5]。通常，地下水资源评价是指对地下水资源允许开采量的评价。1.2地下水模拟的研究现状1.2.1国外研究现状最早的地下水动态分析是基于直观和最简单的水平衡方法，以与水文地质学方法相比较[6]。19世纪中期，亨利·达西在总结前人实践的基础上，提出了水在多孔介质中的线性渗透定律，即达西定律。后来，J.Dupuit在达西定律的基础上，研究了单向和平面径向稳定运动，奠定了地下水稳定流理论的基础。1905年，E.Mailer首次用解析法论证了泉水流量的预测方法[7]。1935年，C.V.泰斯提出了地下水流向承压水井的非稳定流公式，开创了现代水文地质计算的历史。20世纪50年代，随着深层承压水的开发利用，C.E.Jacob、M.S.Hautush等人研究了溢流补给的情况，随后出现了考虑无压含水层后期反应和不完全井的解析解。同时，将在稳定流计算中行之有效的叠加原理和映射方法应用于非稳定流计算中，以解决井组干扰、边界影响、抽水流量阶跃变化等非稳定流问题。20世纪50年代后期，T.H.Kamensky在用解析法分析群孔地下水动态的基础上，系统地研究了降水入渗条件下的有限差分法，并用它来预测地下水动态的变化。20世纪60年代后期，随着计算技术的进步，数值模拟方法应用于地下水计算，地下水的理论和工程分析取得突破性进展。先后出现了二维流动平面(剖面)模型、准三维流动模型、三维流动模型和耦合模型[8-12]。近年来国外出现的几个主要模型，如R.Bravo(美国)的休斯顿模型和美国弗吉尼亚州的滨海平原三维模型，确定了地下水在无约束系统中的运动方式、主要补给区和排泄区以及地下水的循环次数。澳大利亚昆士兰大学系统地研究了地下水动力学的数值解和建立模型的方法步骤，并在实践中模拟了柯利盆地地下水与地表水的动态耦合；模拟了海水入侵引起的滨海含水层污染物迁移过程，并将数值模型应用于含水量计算、开采预测、管理规划和含水层修复，取得了显著效果。目前地下水系统的数值模拟方法[13]主要有有限差分法(FDM)、有限元法(FEM)、边界元法(BEM)和有限分析法(FAM)。国外这方面的研究主要是针对数值模拟方法的弱点，提出新的思维方法，采用新的数学工具，分析不同尺度下的变化，合理描述地下水系统量的不确定性和模糊性因素。按照Anderson等人[14]提出的工作程序，要建立正确的、有意义的地下水系统数值模型，需要做以下工作:确定模型目标、建立水文地质概念模型、建立数学模型、设计求解模型、校正模型、分析校正敏感性、验证预测模型、分析预测敏感性、设计模型并给出模型结果、检查后续重新设计模型。Wood[15]根据Newman猜想提出了二维地下水运动有限元计算的时间步长条件。Kim[16]等人数值模拟了地下水抽取引起的Noordbergum效应(反向水位流态化)现象，并阐述了其机理。Scheibe[17]等人分析了不同尺度下的地下水流及其迁移行为。Ghassemi[18]指出，三维模型可以详细解释含水层系统的三维边界条件和抽水应力情况，而二维模型无法妥善处理。李曙光等人[19]指出数值模型不能解决预测中不确定因素的问题，创造性地提出了随机地下水模型，可以解决小尺度过程的均值分布和不同尺度的问题。Mehl[20]等人提出了二维局部网格细分法的有限差分地下水模型，提供了新的插值和误差分析方法，提高了模拟结果的可靠性。1.2.2中国的研究现状中国的地下水数值模拟始于20世纪70年代末。随着地下水数值模拟研究的发展，地下水数值模拟已经成为研究地下水运动规律和定量评价地下水资源的主要方法，其发展趋势已经远远超出了原来作为计算方法的领域，成为模拟水文地质过程、发生和发展演化规律的主要方法。目前，我国已经建立了国际地下水模拟中心(IGWMC)P.VenderHeijde分类中的所有模型，即预测模型(包括水流模型、热量运移模型、变形模型、多目标模型)、管理模型和识别模型。研究内容涵盖了饱和带、非饱和带和饱和-非饱和带，基本满足了我国国民经济建设和发展的需要。随着非稳定流理论的发展和计算机的广泛应用，各种复杂条件下的地下水运动都可以用数值方法求解。近年来，随着新技术、新方法的广泛应用，我国该领域的科学家做了大量工作，发现了地下水系统数值模拟模型建立中的问题，并通过数值模型理论与相关研究方向的结合，在理论和方法上不断创新，不断提高模拟结果的可靠性。陆文熙[21]对地下水运动数值模拟中的边界条件进行了分析，提出模型预测前应考虑自然因素、人类活动和邻近水流条件的耦合作用，并先预测边界条件。薛玉群等人[22]介绍了Ms-FEM(多尺度有限元法)的基本原理，并将其应用于非均质多孔介质中的流动问题。通过比较计算结果，得出多尺度有限元法比传统有限元法更有效的结论。该校地球科学系[23]用伽辽金有限元法建立了柳林泉平面二维区域地下水渗流模型，描述了柳林泉域的地下水渗流，模拟效果良好。在此基础上，预测了柳林电厂水源对区域地下水流场和柳林泉的影响。工业大学[24]根据市场水源地裂隙岩溶水的赋存和运移特征，建立了三维有限元数学模型，采用等参有限元法求解，并根据大规模群孔实验数据，识别了模型参数，进一步评价了裂隙岩溶水资源开采方案，为水源地开采提供了依据，对类似条件地区的地下水资源评价具有重要的参考价值。中国地质科学院水文地质与环境地质研究所[25]在黄河三角洲利津地区建立了三维地下水数值模型，分析预测了地下水动态场的变化和浅层地下水淡水开采引起咸水入侵的可能性。中国地质大学[26]在分析水文地质条件和建立黄河下游影响带数学模型的基础上，用FEFLOW软件建立了三维地下水流模拟模型，详细描述了模型中防渗墙的处理，并用模型模拟了黄河侧向渗流和防渗墙建成后地下水流场的变化。1.3课题的提出、研究内容、方法和目的。1.3.1提出主题本课题的研究是结合正在施工的长江大桥特大沉井施工降水方案的研究进行的。长江大桥是世界上第一座三塔两跨的悬索桥。长江大桥东距江阴长江大桥57公里，西距润扬大桥66公里。是省“五纵九横五连通”高速公路网的重要组成部分。大桥全长62.088公里，由北接线、跨河主桥、夹江大桥和南接线四部分组成。工程起于童宁高速公路宣堡路口，在永安洲镇过长江，向西在扬中镇夹江泡沙，经桥镇进入境内，止于沪宁高速公路唐庄路口。跨河主桥为主跨1080m×2的三塔双跨钢悬索桥。主桥和夹江大桥总长9.726km，桥面宽度33m。主桥通航净空高度不小于50m，净宽不小于760m，可满足5万吨级巴拿马货轮的通航要求。夹江大桥通航净高和宽度分别不小于18m和100m。1.3.2研究能力南北锚碇是悬索桥的重要组成部分，锚碇结构采用沉井。根据实际施工情况，本文仅对北锚碇施工降水方案进行分析研究。北锚碇设计尺寸分别为52.4米×68.3米，基底高程为-55.0米[27]。锚固的截面图如图1-1所示。图1-1锚固的横截面图根据沉井施工方案，沉井施工过程中需要减少用水量。根据地质资料，锚地被松散沉积物覆盖，沉积物由第四系和上第三系组成，土层富含水分，与长江关系密切。沉井深厚的降水深度会引起沉井周围的地面沉降和开裂，严重威胁周围建筑物的安全和稳定，尤其是附近长江大堤的安全。该堤防为一级堤防，其安全至关重要。因此，有必要对沉井施工中降水后的土体变形、渗流变形和堤防整体稳定性进行分析评价[28]。在地面沉降和堤防变形的预测和控制中，需要综合研究地层的空间分布形态、结构特征、物理力学性质、地下水分布和抽水状况。因此，利用现代计算机信息技术建立地层的可视化模型，直观地再现锚固区不同性质(地质成因、分布方式、物质组成、结构特征、物理力学性质、水力特征)地层的空间分布格局以及地层之间的组合关系，是十分必要的。1.3.3研究方法目前最常用的数值模拟方法是有限元法和有限差分法[13]。根据不同的方法，国外编制了不同的数值模拟软件，主要有加拿大Waterloo水文地质公司的VisualMODFLOW、德国WASY水资源规划与系统研究所的有限元软件FEFLOW、美国杨百翰大学环境模型系统实验室的GMS等。[29]这些国际流行的地下水流模拟软件有其共同的特点，即:功能齐全的数据输入和输出手段；友好的人机界面，可视化效果好。本文利用有限元软件GMS(地下水模拟系统)对长江公路大桥沉井施工的地下水流路径进行三维模拟计算，分析河堤的渗透稳定性以及河堤和地面的沉降变形。GMS[30]是美国杨百翰大学环境模型研究实验室和美国陆军排水工程试验工作站在整合MODFLOW、FEMWATER、MT3DMS、RT3D、SEAM3D、MODPATH、SEEP2D、NUFT、UTCHEM等现有地下水模型的基础上开发的地下水模拟综合图形界面软件。其图形界面由下拉菜单、编辑栏、常用模块、工具栏、快捷键和帮助栏六部分组成，使用起来非常方便。目前国内利用该软件分析地下水运行系统并在实际工程中实施的还不多见。下面简单介绍一下本文所需的GMS模块。MODFLOW[31]是美国地质调查局于20世纪80年代开发的一套三维有限差分数值模拟软件，专门用于多孔介质中的地下水流动。自MODFLOW问世以来，因其具有程序结构模块化、离散方法简化、求解方法多样化等优点，被广泛用于模拟井流、河流、泄流、蒸发和补给对非均质复杂边界条件水流系统的影响[32-34]。MODPATH是一个三维示踪模型，用于确定给定时间内稳定或不稳定流中颗粒的迁移路径。它与MODFLOW一起使用。根据MODFLOW计算的流场，MODPATH可以追踪一系列虚拟粒子，模拟污染物从用户指定的地方溢出的运动。这种追踪方法可以用来描述一口井在给定时间的拦截面积。SEEP2D是用于计算坝段渗流的二维有限元稳定流模型。它可用于模拟承压和自由表面流动问题，也可用于模拟饱和区和非饱和区的水流。对于自由表面流动问题，模型只能限于饱和带。根据SEEP2D的结果，可以制作一个完整的流网。PEST和UCODE是用于自动参数调整的两个模块。自动估计参数时，交替使用PEST或UCODE调整所选参数，并重复用于MODFLOW、FEMWATER等的计算。直到计算结果与现场观测值一致。MAP使用户能够快速构建概念模型。在地图模块下，以TIFF、JEPG等地图为底图，在地图上确定代表源、汇项、边界、含水层不同参数区域的点、曲线、多边形的空间位置，快速建立概念模型。tin(不规则三角网)通常用于表示相邻地层的界面，多个tin可用于建立实体模型或三维网格。实体是不规则三角网(TIN)建立后，通过一系列操作产生的实际地层的三维模型。钻孔数据(钻孔数据)用于管理两种格式的钻孔数据:样本和地层。样本数据用于制作等值面和等值线，地层数据用于构建TIN、实体和三维有限元网格。1.3.4技术路线地下水数值模型的建立需要考虑地下水含水层系统中含水层、弱透水层和隔水层的空间组合分布，以及水文地质参数的时空变异性。根据本次研究的目标和要求，制定了如图1-2所示的技术路线。图1-2论文技术路线1.3.5达到的目的此次通过模拟建立三维可视化地层，动态直观地展示水流的运动状态和路径，为降水施工方案提供依据；地层的不均匀性导致了地面的不均匀沉降。在模型分析过程中，可以预测周边地区和河堤的沉降，有利于控制地面沉降的发生和发展。河堤的渗流模拟将显示河堤断面完整的流网和等势线，可以准确分析河堤的渗流稳定性。1.4本章概述本章主要介绍了地下水模拟的重要性和国外的研究现状，并根据研究区的实际情况提出了论文的研究方法、技术路线、目的和创新点。3地下水系统的数值模拟3.1地下水系统概化模型要模拟地下水流，首先要确定研究区域。研究区域是客观存在的，客观存在本质上是一个连续的整体。只有准确了解研究区的成因(地质构造和水文地质条件)和外因(边界条件)，才能将其视为模拟研究的一个整体。数值模拟的目的是再现客观存在，研究其客观规律[36]。但实际水文地质条件非常复杂，很难建立数值模型来描述计算区的地下水系统。因此，实际水文地质条件应根据水文地质条件和地下水资源评价目的进行简化。这个过程称为水文地质条件概化，其原则是:根据评价的目的和要求，概化的水文地质概念模型应反映地下水系统的主要特征；概念要简单明了；该概念应能用于进一步的定量描述，从而建立符合研究区地下水运动规律的微分方程的定解。水文地质条件概化通常包括以下几个方面:计算区几何形状概化；含水性质的概括；边界性质的推广；参数性质的概括(同质或非同质、各向同性或各向异性)；地下水流动状态的概化(一维、二维或三维)[37]。3.1.1地下水系统的含水层结构特征[35](1)第一含水层(潜水)属全新统含水层，岩性为灰黄色，底部为灰壤土、淤泥、淤泥质壤土，厚度6.6~8.3m，透水性一般，富水性较差，主要接受大气降水和场区沟水补给。根据抽水资料，抽水过程中下层水位无明显变化，说明第一层地下水与下层地下水无明显水力联系，溢流影响不明显。(2)第二含水层(承压水)该含水层属全新统，主要岩性为粉砂，局部为细砂。顶部标高为-4.01~-5.76米，底部标高为-23.21~-23.66米，层厚17.90米~19.20米，砂质纯净，透水性好，富含水。通过两组抽水数据的对比可以看出，靠近长江的试验井涌水量较大。(3)第三含水层(相对较弱的含水层)该含水层属于全新统，主要岩性为亚粘土夹粉砂或粉砂，局部为粉砂。顶部标高-23.21~-23.66m，底部标高-49.46~-51.91m，层厚27.80m~28.70m，单井涌水量较小，受长江潮汐影响不明显。抽水时，降落漏斗的影响半径较小，观测井的水位落差较小，因此该含水层的透水性较弱。3.1.2地下水系统的水交换特征从空间上看，地下水流整体以水平运动为主，垂直运动为辅，地下水系统符合质量和能量守恒定律。含水层分布广、厚度大，常温常压下地下水运动符合达西定律。考虑到浅、深层之间的流动交换和软件的特点，地下水运动可以概化为三维流动。地下水系统的垂向运动主要是层间溢流，三维地下水模型可以很好地解决溢流问题。地下水系统的输入和输出是随时间和空间变化的，所以地下水是不稳定流。参数随空间变化，反映了系统的异质性。根据抽水试验资料，含水层具有明显的方向性，因此含水层可以概括为非均质性和各向异性。因此，根据研究区水文地质条件和地下水补给及动态变化特征，将研究区概化为一个非均质、各向异性、三维空间结构和不稳定的地下水流系统。3.1.3地下水系统的边界条件边界条件的研究是建立水文地质三维概念模型的一个重要方面。边界可分为外边界和外边界。外边界是指每个含水层组的周围边界，边界是指连接多个含水层的研究区域的边界。根据本区的水文地质条件，研究基本上构成了一个完整的地下水补给、直径和排泄系统。主要有两个边界:一是侧边界，南部沿长江，接受长江强大的供水。在模型中，可以设置为已知的水头边界，水位值为长江边界附近的水位；东、北侧地下水等水位线基本垂直于边界，故东、北侧简化为零通量边界(流线边界)。第二个是垂直边界。潜水含水层的自由水面是该系统的上边界。潜水通过这一边界与系统外部进行垂直交换，如长江补给、河塘补给、大气降水入渗补给、蒸发排泄等。由于研究区建设时间短，不考虑大气降水入渗补给和蒸发排放。3.2地下水系统的数学模型数学模型的表达确定一个水文地质问题的数学模型，需要从概念模型出发，用一组数学关系来描述其数量关系和空间形式，从而反映所研究的地质体的地质和水文地质条件以及所要描述的地下水流的基本特征，达到再现或复制一个实际地下水系统的基本状态的目的。这组数学关系通常被称为数学模型[38][39]。这组数学关系式不仅概括了已知部分(水头等。)的系统，还包含未知部分(参数)。通过已知数据，即数学模型的解，预测系统的性质。在大多数情况下，地下水的动态迁移是用三维偏微分方程来描述的。这些方程，加上附加的边界条件和初始条件，构成了一个数学模型[36]。对于承压水和非承压水的流动运动，在不考虑水密度变化的情况下，以水头为基本变量，多孔介质中地下水流动的有限差分公式数学模型可以用下面的偏微分方程[38][40][41]表示:（3-1）（3-2）（3-3）（3-4）其中：——渗透系数在方向上的分量，我们假定渗透系数的主轴方向与坐标轴的方向一致；——水头；——单位体积流量，代表流进汇或来自源的水源；——孔隙介质的贮水率；——时间；——研究区各层初始水头值；——研究区各层第一类边界上的实测水头值；——研究区各层第二类边界上的单位面积流量；——一类与二类边界；——承压含水层边界上的法线。公式(3-1)加上相应的初始和边界条件构成了描述地下水流系统的数学模型。该数学模型的解是描述水头值分布的代数表达式。在规定的空间和时间内，求得的水头H应满足边界和初始条件。数学模型的求解数学模型建立后，其求解包括正演模拟和反演。如果给定含水层的水文地质参数、边界条件和初始条件，就可以求解含水层在各个时刻的水头值，预测水头。这类问题称为正演问题，即正演建模。如果水文地质参数是由动态观测资料或抽水试验资料确定的，这个过程称为反演，即前者的反问题或参数反问题。为了验证所建立的数学模型的正确性，求解中一般需要进行模型辨识，即反演过程。模型经反演验证后，利用反演得到的水文地质参数进行正演模拟，求解所需项目，如水位预报、涌水量预报等。上述数学模型很难用解析方法求解，但通常借助计算机用数值方法求解。本文采用GMS软件进行数值模拟，其数值方法为有限差分法。有限差分法[42][43]的基本思想是将研究区域用差分网格离散化，对每个节点用差商近似控制方程中的微信业务，然后耦合初始条件和边界条件，求解封闭的线性代数方程组。在用有限差分法求解的过程中，连续的时间和空间被分成一系列离散的点。在这些点上，连续偏导数也被水头差公式所代替。所有的未知点组合起来，这些有限差分公式就形成了一个线性方程组。然后联立求解线性方程组，使得得到的解是各离散点水头的近似解。数值解虽然不能给出一个描述水头随时间和空间变化的代数表达式，但可以用来解决很多实际问题。3.3地下水系统的数值模型模型的预处理3.3.1.1电网公司在数值模拟中，研究区域将由节点和相关的有限差分网格或有限元组成的离散域来代替。节点间距的选择是网格设计的关键步骤。在水头变化较大或需详细了解水头分布的地方，节点间距应较小，但相邻节点间距或机组大小不应相差太大，应逐渐过渡。观测井和抽水井尽量放在节点上，这样提供的水头和流量更准确，也更方便校核[38]。本研究区计算面积为181350，采用GMS软件进行自动三角形网格剖分，剖分单元13058，结点9338。由于在河流和抽水井处水头变化比较大，因此对网格进行了加密。网格划分见图3-3。3.3.1.2时间离散时间步长的选择和空间步长的选择一样是模型设计的一个关键步骤。在抽水试验时，水头下降速度不均匀，经历着一个由迅速下降至逐渐缓慢的过程，因此，抽水初期用较小的时间步长，以后随着水头下降速度的减缓，逐步加大时间步长。时间步长的增大一般按下式进行：，deMarsily[44]认为最好的选择是应用,本文所用的时间步长的增大系数为。3.3.1.3水文地质参数及初始值的确定水文地质参数是表征含水介质蓄水能力、放水能力和地下水运动能力的指标。水文地质参数的分析、研究和计算是地下水资源计算和评价的重要环节之一。参数选择的准确性直接影响地下水资源计算和评价的准确性。水文地质参数往往是根据地下水动态观测资料，通过野外试验、室内试验及相关理论公式计算得出，或者采用数值反演的方法获取参数。本次模拟采用GMS中的LPF（LayerProperyFlow）子程序包，需要输入每一含水层的水平向渗透系数和垂直向渗透系数。在非稳定流计算中还需要输入贮水率和给水度。[41][45][46]水平和垂直渗透系数的初步确定采用了研究区现有钻孔地质资料和抽水试验资料获得的实测参数值。然后基于这些离散点的值，通过特定方法的空间区划方法形成研究区的初始参数值，通过模型调试确定模型的最终参数值。（1）贮水率弹性贮存或弹性释放是含水层的重要水文地质特征。为了从数量上评价，需要表明含水层中贮存水量的变化和承压水头相应变化之间的关系，通常采用贮水率或贮水系数来表明这种关系。贮水率表示当含水层水头降低一个单位时，由于含水层骨架的压缩和水的膨胀而从单位体积含水层中释放出来的总水量。当承压含水层的厚度固定不变时，有，它所释放的水量完全是弹性水量。潜水贮水系数所释放的水量包括两部分：一部分是含水层由于压力变化所释放的弹性水量；另一部分是水头变化一个单位时所疏干含水层的重力水量，这一部分水量正好等于含水层的给水度，由于潜水含水层的弹性变形很小，可以用给水度近似代替贮水系数。贮水系数是没有量纲的参数，其确定方法是根据非稳定流抽水试验资料，利用配线法（标准曲线对比法）的降深—时间配线法推求[41]。由式：（3-5）（3-6）式中：——抽水井流量；——抽水井井壁降深；——抽水井半径；——抽水时间；——含水层的贮水系数；——导水系数；——井函数；对式（3-5）、（3-6）取对数，得（3-7）（3-8）从解析几何中知道，如果有，，那么曲线和的形状是一样的，只是曲线相对于曲线在横坐标上位移了，在纵坐标上位移了。因此，如果在双对数纸上做的关系曲线，在另一模数一样的双对数纸上做关系曲线，显然，这两条曲线的形状也是一样的，只是坐标位移了一个距离而已。降深—时间配线法的方法如下：根据数值表，在双对数纸上作一的关系曲线，称之为标准曲线。在另一模数一样的透明双对数纸上，根据抽水资料作关系曲线（即降深—时间关系曲线），称之为实际资料曲线。把实际资料曲线重叠在标准曲线上，保持坐标轴平行，移动到二根曲线完全重合为止。在重合的双对数纸上任选一点作为配合点，读出改点在两对数纸上相应的坐标，，和的值，再把值代入公式（3-5）、（3-6）中，即可解出和来。具体计算见表3-1。表3-1贮水率计算表井号配比点(h)(m)S1-110.08331.358.11.776.25.082.93×10-31.577×10-420.2782.53034.3632.1113.920055.30S1-210.08330.281.50.42.8820.2220.7412.8831.6672292.93.02S4-110.20.44100224.761.3820.60.930041.35341.1200051.38S4-210.11670.0830.42.5920.19440.21512.4731.390.423022.41S5-110.10.24100062.940.0954.1×10-41.451×10-5220.32100040.953100.56500070.95S5-2110.2850050.158220.29100050.1523100.34500060.156（2）给水度给水度是表征潜水含水层给水能力或储水能力的一个指标[46][47]。它与岩性有关，随排水时间、潜水埋深、水位变化幅度与水质的变化而变化。本次模拟的给水度根据研究区含水层岩性与颗粒组成，取给水度为0.06。(3)初始条件初始条件是在选定的初始时间渗流区的水头分布。初始时刻一般理解为根据需要任意选择的时刻，但不一定是泵送的实际开始时刻。本文利用各观测井此时的实测水头，通过插值得到模拟初始水头。在动态稳定流条件下，流入系统的水等于流出系统的水，水头在空间上是变化的。Anderson和Woessner认为“选择从已识别模型获得的稳态水头解作为初始条件是一种标准做法”。[48]因此，在具体模拟时，选取一段相对稳定的水位作为初始时间，输入边界条件，得到模型的稳态解，可作为下一次非恒定流模拟的初始条件。建立拟合期初始流场时，采用2006年12月13日8:00至2006年12月17日20:00的实测水位，通过线性插值得到同一水位上各单元中心点的初始水位值，然后将模拟的初始水位值赋给模型。在模型预测中，选取一段相对稳定的水位作为初始时间，输入边界条件，将建立的模型作为下一次非恒定流模拟的初始条件，得到短时间内的稳态解。在3.3.1.4加工源和汇项目1.增刊研究区地下水的侧向补给主要来自河流渗漏补给和研究区周边边界的补给，垂向补给主要是含水层溢流补给，调用GMS中的河流子程序包。2.排泄项目地下水的排泄项主要以抽水井的形式表示，调用GMS中的井的子程序包，在模型中用节点表示，并给定一定的流量，负值表示抽水。3.3.2建立三维可视化地层模型信息系统。根据研究区含水层的埋深、厚度和岩性，构建了GMS地下水数值模拟模型。建模过程如下:(1)根据地质钻孔数据，由钻孔坐标确定计算区域内各钻孔的坐标位置，生成钻孔柱状图。(2)将地面和河堤的实际位置和尺寸导入到模型中。(3)根据TINS模块，在井眼内用反距离加权差分法生成固体地层，如图3-1和3-2所示。(4)利用网格法将实体模型转化为网格模型。这一步实现了有限差分方程中计算系数的转换。网格是指地下水的建模过程和模型识别验证都是直接在三维网格上操作的。网格法是最简单、最基本的建模方法，对于水文地质条件简单、研究面积小的模型是一种有效的方法。网格建模是基于对网格的直接操作，因此参数的赋值、边界条件的指定和参数的调整都非常繁琐，工作量相当大，枯燥且容易出错。在这个建模过程中，采用网格方法和概念模型联合建模。概念建模的思想是用覆盖层中对应实体的属性来表示三维网格的属性，通过模型赋值和参数调整来直接操作覆盖层上的实体。一个覆盖可以包含点(点，代表曲线的顶点；节点，代表曲线的终点/节点)、圆弧(代表河流、沟渠等线性实体)、多边形(代表渗透系数分区等平面实体)。这些实体有独特的识别号码，可以分配属性[49]。该模型的计算网格模型如图3-3所示。(5)在三维网格下设置Modflow的初始条件，给定井的抽水量、河流和边界水位的变化，选择计算包和迭代计算方法。在MAP模块下，将MAP文件中的数据转换成Modflow，实现了有限差分方程中源汇项的转换。(6)检查建立的模型和给定的求解条件是否正确，不断调整修正，然后运行Modflow进行计算。图3-1研究区域的物理模型图3-2研究区固体地层剖面图图3-3研究区网格模型3.3.3型号识别和检验通过对研究区水文地质的合理概化，建立了水文地质概念模型，并根据这一物理概念模型建立了相应的地下水渗流数学模型。经过一系列的高度概化和相应的模型转换，虽然每个阶段都有理论依据，但转换后的最终模型能否全面客观地代表研究区具体的水文地质条件和特征？这需要模型识别和验证。模型识别和验证的过程是整个仿真中极其重要的一步[50]。模型的鉴定是为了证明所建立的模型确实具有重现现场实测水头和流场的能力，并且模型所反映的这些场是符合实际情况的[51]。因此，需要对模型进行检验，即将模型的预测结果与实际观测曲线或通过某种实验或对含水层的某种影响而得到的某一地区地下水动态长期观测数据进行比较，看是否一致。如果不是，模型应该被修正。通常通过反复修改参数，调整一些源项和汇项，就可以达到理想的拟合效果。模型经过识别和校正后，才能初步说它真正代表了所研究的地质体，是实际水流系统的复制品，才能用于进一步的计算或预测。辨识模型参数有两种基本方法:(1)试错法，手动调整参数；(2)自动调整参数。该模型的识别和检验过程采用试评价和修正的方法[52][53]。该方法预先给定一组参数的估计值，输入模型进行计算，比较计算结果与实测结果的误差。如果误差不满足精度要求，则调整刚输入的一组参数，直到达到一定精度，停止计算。此时，使用的一组参数值将被视为实际参数值。本文数值模拟结果的可靠性是计算地下水位与实测地下水位的拟合程度。模型的识别和验证主要遵循以下原则[54-56]:①模拟的地下水流场应与实际的地下水流场基本一致，即模拟的地下水等高线应与实测的地下水位等高线相似；②模拟地下水的动态过程应与实测基本相似，即模拟与实际地下水位过程线形状相似；③从平衡的角度来看，模拟的地下水平衡变化应与实际情况基本一致；④确定的水文地质参数应符合实际水文地质条件。根据上述四个原则，对研究区的地下水系统进行了识别和验证。通过反复调整参数，识别水文地质条件，确定模型结构和参数。本研究由于缺乏足够的可以控制整个区域的长期水位观测数据，无法通过地下水的长期动态数据进行模型识别，因此采用现场抽水试验期间的水位动态观测数据，并将识别期定为抽水试验期。将抽水试验中各观测井的水位观测数据代入模型，观测各水位观测井的计算水位与观测水位的拟合，作为调整参数的依据。GMS提供了直观的校准标记功能，通过该功能可以在模型的平面图中直接查看拟合效果，如图3-4所示。符号定义如下:当模型在观测井上得到的误差值(计算值与观测值之差)在设定的允许范围内时，符号显示为一个短的绿色条。随着误差的增大，它的颜色会逐渐由绿变黄变红，长度也会变长。标记还可以给每个观察井赋予权重，以表示每个观察井对拟合效果的贡献不同。利用刻度线，可以直观地看到各观测井水位拟合的总体情况，便于参数调整。图3-4校正标记刻度线只能用来大致了解研究区域的水位拟合情况。要想定量了解拟合效果，就得计算观察水位对比图。本文以各含水层的水平渗透系数和水平渗透系数与垂直渗透系数之比作为待调整参数。在模型辨识过程中，根据设定的允许误差和模型辨识原理，根据实测水位值反复调整辨识不同的参数，并给出两个参数的周边值。水位识别结果见图3-5~图3-12。从图中可以看出，由辨识出的参数计算出的模拟水位线与实测水位线趋势相同，误差在允许范围内。图3-5第一含水层1号观测井图3-6第一含水层2号观测井图3-7第二含水层1号观测井图3-8第二含水层2号观测井图3-9第二含水层3号观测井图3-10第二含水层4号观测井图3-11第三含水层1号观测井图3-12第三含水层2号观测井3.3.4模型识别和评估模型辨识通过比较实测水头和模拟水头等值线图，提供了一种可视化的定量方法，给出了辨识误差空间分布的一般信息。而根据野外数据绘制的等值线图含有绘图引起的误差，不适合作为识别的唯一证据，需要对识别结果进行定性和定量评价。测量头部和模拟头部的散点图是显示拟合程度的另一种方法。每一点对直线的偏离应该是随机分布的。用列表法表示水头观测值、模拟值以与它们间的误差值和某种误差值的平均是展示识别结果的另一种常用方法。误差值的平均就作为识别过程中的平均误差。识别的目的就是使这个误差最小。一般采用均方根（）误差[38]来表示水头观测值和模拟值间误差平方和的均值：式中：为识别值的数目。根据模型识别的结果，对各参数组合进行识别评估，各含水层的见表3-2～表3-4。从表中所列的，选择最小误差的参数组合作为各含水层的参数值，具体值见表3-5。表3-2潜水含水层模型参数识别与评价汇总表Kh=1.0，Kh/千伏=10Kh=0.9，Kh/千伏=8Kh=0.8，Kh/千伏=6Kh=0.8，Kh/千伏=7Kh=0.7，Kh/千伏=5Kh=0.7，Kh/千伏=6Kh=0.6，Kh/千伏=4Kh=0.5，Kh/千伏=3S3一号0.1080.10430.10630.0910.09610.0770.0880.0851S3-20.0420.04460.05040.0420.05040.04040.0530.0594总数0.150.14890.15670.1330.14650.11740.1410.1445表3-3相对弱含水层模型参数识别与评价汇总表Kh=6.8，Kh/千伏=24Kh=6.9，Kh/千伏=24Kh=7.0，Kh/千伏=24Kh=7.1，Kh/千伏=24Kh=6.8，Kh/千伏=22Kh=6.9，Kh/千伏=22Kh=6.7，Kh/千伏=20Kh=6.8，Kh/千伏=20Kh=6.6，Kh/千伏=18S4一号0.31780.28110.24080.20740.28460.2490.29330.25690.2983S4-20.33620.35070.37550.38990.37150.3860.38340.39750.3987总数0.6540.63180.61630.59730.65610.6350.67670.65440.697表3-4承压含水层模型参数识别与评价汇总表Kh=21.4，Kh/千伏=6Kh=21.6，Kh/千伏=7Kh=21.6，Kh/千伏=10Kh=21.8，Kh/千伏=8Kh=21.8，Kh/千伏=10Kh=21.8，Kh/千伏=11Kh=22，Kh/千伏=8Kh=22，Kh/千伏=10Kh=22，Kh/千伏=12Kh=22.2，Kh/千伏=16Kh=22.2，Kh/千伏=18Kh=22.4，Kh/千伏=22Kh=22.6，Kh/千伏=22S1一号0.25360.18570.21120.12610.13960.14660.08060.08660.09380.08230.08310.09560.1496S1-20.03740.07620.060.12330.11050.10580.17580.16340.1530.19180.18190.21750.2668总数0.2910.26190.27120.24940.25010.25240.25640.250.24680.27410.2650.31310.4164Kh=23，Kh/千伏=8Kh=23.2，Kh/千伏=8Kh=23.2，Kh/千伏=10Kh=23.4，Kh/千伏=8Kh=23.4，Kh/千伏=14Kh=23.6，Kh/千伏=8Kh=23.6，Kh/千伏=10Kh=23.6，Kh/千伏=14Kh=23.8，Kh/千伏=10Kh=23.8，Kh/千伏=14Kh=23.8，Kh/千伏=20Kh=24，Kh/千伏=14Kh=24，Kh/千伏=18Kh=24.2，Kh/千伏=16Kh=24.2，Kh/千伏=18S4一号0.15260.1530.15030.170.150.19970.18550.16370.21950.1910.16160.22630.19920.2510.2364S4-20.140.1130.1240.10050.1240.10.09850.1060.1020.09780.1140.1010.09740.10680.1014总数0.29260.2660.27430.27050.2740.29970.2840.26970.32150.28880.27560.32730.29660.35780.3378表3-5含水层参数含水层名称（）潜水含水层0.76-0.06承压含水层22.581.577×10-4-相对弱含水层7.0241.451×10-5-3.3.5模型检验由于辨识中的不确定性，由辨识模型得到的一组参数可能不能准确地代表现场的实际值。为了进一步验证所建立的数学模型和模型识别后确定的水文地质参数的可靠性，有必要利用现有的另一时期的地下水位动态观测数据对数学模型进行检验。模拟值与现场实测值应在预设的允许误差范围内保持一致，且允许误差应小于拟合计算时水位变化值的10%。当水位变化值较小时(小于5m)，水位拟合误差一般应小于0.5m[38]如果所有结果显示一致，则模型可靠、实用，可用于预报。如果不是，必须对辨识得到的最后一组参数进行修正，重新进行辨识和检验，直到新的一组参数显示模拟值和测量值都在预设的允许误差范围内一致。由于缺乏该研究区域的长期观测数据，试验由另外三组独立的稳定流数据完成。图3-13第一含水层计算水位与实测水位过程线对比图3-14第二含水层1的计算水位与实测水位过程线对比图3-15第二含水层2的计算水位和实测水位过程线对比图3-16第三含水层计算水位与实测水位过程线对比通过对计算地下水位(流场)与实测地下水位(流场)拟合误差的统计，发现在模型的识别检验阶段，所有观测井水位拟合良好，模型计算值与观测井反映的实测值变化趋势基本一致。拟合误差小于0.5m，能较好地反映该点的水位动态趋势，如图3-13~图3-16所示。结果表明，含水层结构、边界条件的概化、水文地质参数的选取以及源汇项的确定都是合理的。所建立的数学模型能够真实地描述研究区地下水系统的特征，具有较强的模拟性。建立的数值模型可用于地下水资源评价。经分析，误差的主要来源是各种源、汇项的统计误差，野外观测的水文地质参数的误差，以及水文地质边界的确定，往往具有较大的人为因素和不确定因素。另外，地层数据的准确性导致了地层模拟的误差。还需要注意的是，模型算法采用迭代求解，迭代法得到的解只是差分方程的近似解。它的精度受许多因素的影响，如选定的收敛指数和迭代方法本身。即使每个时间段内得到的解是精确的，也只是相对于那个时间段内建立的差分方程而言。对于偏微分方程公式，其数值解只是一个近似解。因为有限差分法得到的解与相应的解析解相比总是存在截断误差。一般来说，这种截断误差会随着网格间距的增大而增大。即使得到了基本偏微分方程的精确解，对于现场条件来说，它仍然是一个近似解。模型预测通过对模型的识别和检验，表明所建立的模型可以作为真实地下水系统的代表，可以用来预测未来的地下水运动状态。所谓预测，就是将规划的地下水开采量和预测的地下水源和汇量代入已识别和检验的数学模型，通过模型的运算计算，计算出不同降水深度下沉井内地下水位和流场的响应。3.3.6.1一类边界的预测一级边界的水位受自然因素控制。在预测阶段，在模型预测之前预测边界条件。由于桥址无水文站，缺乏水文资料，根据相邻站的水文资料计算多年平均水位和五年水位，分别为3.0m和5.5m。由于大桥位于感潮河段，以上水位均为高潮位，中间水位的计算由程序自行插值完成。3.3.6.2模式中源汇项的预测在进行模型预测之前，需要对模型中的源项和汇项进行预测，特别是抽水量。根据沉井施工方案，抽水井等间距分布在沉井周围，每侧5口，共16口。抽水井与沉井之间的距离为10m。抽水井布置见图3-17。图3-17研究区施工降水井布置图3.3.6.3预报周期的确定根据沉井施工方案，施工降水分为四期，每期持续一个月。预测时段见表3-6。表3-6预测期划分预报期ⅠⅡⅢⅣ经历的时间2008-3-1~2008-4-12008-4-1~2008-5-12008-5-1~2008-6-12008-6-1~2008-6-303.3.6.4模型的预测结果根据施工方案，四期沉井中心降水深度分别为8m、13m、18m和23m。通过模型计算和预测，各抽水井涌水量见表3-7。表3-7机井涌水量(单位:)抽油井标签第一阶段涌水量第二阶段涌水量第三期涌水量第四期涌水量一个2504506808502250450680850三250450680850四250450680850五250450680850六250450680850七170450680850八170230350720九17023035050010170230350500111702303505001217023035050013170230350500141702303507201517045068085016250450680850总进水量32805660857011590根据《南北锚碇水文地质勘察报告》的北锚碇基坑涌水量的估算值，当基坑水位降至-15m时，基坑水主要来源于第二层含水层，基坑涌水量由公式（3-9）计算得：（3-9）式中：——基坑涌水量，；——基坑水位降深，；——潜水含水层厚度，；——渗透系数，——基坑影响半径，；——基坑等效半径，。当基坑水位降至-23m时，基坑水主要来源于第二层含水层以与第三层上层，基坑涌水量由公式（3-9）计算得：将模型预报的涌水量与水位降深作拟合，得出当基坑水位降至-15m时，涌水量，与理论公式计算结果的误差为3.66%；当基坑水位降至-23m时，两种涌水量计算结果的误差为20.55%。该模型考虑了地层分布的非均质性和不同方向渗透系数的差异，考虑的边界条件更接近实际情况。数值模拟预测的基坑涌水量与一般理论公式计算的涌水量一致，误差在可接受的范围内。这进一步验证了模型的可靠性和准确性，为施工的顺利进行提供了可靠的依据。降水过程中，长江附近的水头剖面、上游方向的水头剖面、平面内的水头剖面见图3-18~图3-29。图3-18降水量为8m时长江附近的水头剖面图图3-19降水量13m时长江附近的水头剖面图图3-20当降水量为18米时长江附近的水头剖面图图3-21降水量23m时长江附近的水头剖面图图3-22降雨量为8m时上游方向的水头剖面图图3-23当降水量为13m时，上游方向的水头剖面图图3-24当降雨量为18m时，上游方向的水头剖面图图3-25当降水量为23m时，上游方向的水头剖面图图3-26降雨量为8m时的平面等水位线分布图3-27降水量13m时平面等水位分布图图3-28降水量为18m时的平面等水头分布图3-29降水量为23m时的平面等水头分布图长江纵剖面和平行剖面等值线图中的梯形和粉红色区域表示网格排水的深度和周长。从图中可以看出，随着水位的降低，网格排水的深度和周长逐渐增加，等水位线的变化随着距抽水井距离的增加而减小。形成的降落漏斗水力梯度由外向内由小变大，抽水井区附近水力梯度最大，符合渗流规律。从降水过程线的变化趋势可以看出，靠近长江的沉井两侧水位下降缓慢，周长较小。这是由于扬子江在抽水过程中不断补给，而其他两岸水量不足，所以水位下降快，周长大。3.4本章概述根据本章的水文地质条件和地下水流动态特征，将研究区概化为一个非均质、各向异性、三维空间结构和不稳定的地下水流系统。在此基础上，建立了相应的数学模型。通过钻孔岩性概化，建立研究区水文地质结构模型，从平面和垂向上进行识别，建立三维地质结构模型，结合概念模型建立三维地下水数值模型。然后根据流场对模型进行识别和验证，将识别出的模型应用于地下水流场的预测，预测结果合理可靠。4河堤渗漏4.1渗流理论概述由于长江大桥北锚碇紧邻长江大堤，沉井边缘距江岸仅120m，沉井施工降水将改变江岸的环境地质条件，包括土的物理力学性质、静水压力和渗透力的变化，可能导致滑坡、大变形甚至失稳。一旦滑坡过度变形或失稳，将直接威胁堤防安全，造成严重的经济和社会影响。因此，有必要研究现状及水位下降条件下河堤的稳定性，为岸坡稳定性评价和滑坡防治工程设计提供依据。在影响河堤安全的因素中，渗流起着主导作用。通过渗流计算，可以确定浸润线的位置、动水压力、水力梯度和渗流流量，从而分析堤防在局部渗流情况下的整体稳定性和堤防的安全性。管涌和流土是土的两种基本渗透变形类型。管涌是指土中的细粒土在渗流作用下，在粗粒土形成的孔洞中运动并被带出的现象。管涌可分为深管涌和浅管涌，其中浅管涌的危险性更大。二维地层广泛分布于长江两岸。在这种条件下，上部弱透水层的透水性越低，下部透水砂层中的渗透水压力越高，发生管涌的可能性越大，可能发生在距堤脚数百米的地方。土流是指在向上渗流的作用下，局部土壤表面隆起，或土壤颗粒同时启动并失去的现象。主要发生在基础或下游土坝的渗漏处。在地下水位下降引起的河堤渗漏问题中，渗透力是向下的，所以不会被流土破坏。管涌破坏的形成主要取决于两个基本条件。一方面，土料的颗粒级配，级配不良，特别是缺少中间级配的土，容易发生管涌破坏；另一方面，工程的水力条件是土体管涌破坏的重要原因，在渗流出口处形成较大的水力梯度可能诱发土体管涌破坏。当下游坡脚的水平水力梯度大于堤料管涌的临界水力梯度时，管涌通道将从该点形成，渗流场将发生变化，从而扩大管涌通道周长。可见，地下水位下降过程中形成的动态渗流场削弱了河堤的稳定性，是河堤可能失稳的重要原因。4.2边坡稳定性分析方法概述路堤边坡稳定性分析是工程界普遍关注的问题之一。因为一旦堤防不稳定，就会对两岸人民的生命财产安全造成严重威胁。由此，我国每年汛期都要投入大量的人力物力在河两岸确保堤防安全。水位变化与边坡稳定性密切相关。随着水位的降低，土路堤内形成的渗流场会相应发生变化，从而影响土路堤应力场的变化，产生不利于边坡稳定的不稳定渗流，对路堤影响较大。因此，有必要在水位下降时检查岸坡的稳定性。边坡稳定性分析方法主要有极限平衡理论[57]、塑性极限分析[58]、模糊极值理论[59]和有限元法[60]，大致可归纳为滑动面法和单位应力法两大类。滑动面法主要是以瑞典法、Bishop(1955)法[61]和Janbu，1973)法[62]等为代表的条分法。它可用于均质或非均质土坡的计算，在工程实践中得到广泛应用。单位应力法是应用弹塑性理论计算路堤土体各点的应力分布，然后将某一单位面积的剪应力与其剪应力强度进行比较，来判断其是否安全。这种方法逐渐被采用。此外，结合有限元法，用三角形单元代替垂直条分，能更好地适应各种复杂土层中路堤断面的稳定性分析。终极平衡理论极限方法基于众所周知的摩尔-库仑强度准则，其表达式为式中:为破坏面上的剪应力，kPa；——土的有效粘聚力，kPa；——破坏面上的总应力和有效应力，kPa；——土壤的有效摩擦角，度。极限平衡法的基本原理是:如果土质边坡的稳定安全系数为，当土的抗剪参数(摩擦系数和粘聚力)按倍数减小时，土质边坡最危险滑动面之一的滑动体将处于失稳的极限平衡状态。换句话说，为了得到土质边坡抗滑稳定的安全系数，可以先假定安全系数的取值，用土的摩擦系数和粘聚力除以这个安全系数，作为计算参数进行计算。如果能满足极限平衡条件，假定的安全系数就是我们想要的。否则，重置安全系数并重新计算，直到满足极限平衡条件[63]。极限法的基本特点是:满足力和力矩平衡、莫尔-库仑破坏准则和应力边界条件，通过分析破坏时刻土体的力平衡来解决问题。不考虑土体本身的应力应变关系。在大多数情况下，这个问题是静不定的。处理这个问题的极限方法是引入一些简化的假设，使问题静定可解。这种处理影响了方法的严密性，但对计算结果的精度影响不大，优点是大大简化了分析计算工作，因此在工程中得到了广泛应用。极限分析法1977年，贾政提出了土坡极限分析的两个基本原理，即最大值原理和最小值原理，这是极限分析的理论基础[64]。这一原则认为:(1)如果滑坡可以沿多个滑动面滑动，那么当滑坡失稳时，将沿阻力最小的滑动面破坏(最小原理)。(2)滑坡滑动面确定后，滑动面上的反作用力(以滑坡的力)可以自行调整，从而发挥最大的抗滑能力(最大原理)。同时，专家也指出，这两个原理是相辅相成的，可以成为根据极限平衡原理分析滑坡问题的有力武器。极限分析法采用刚塑性假设，其应力应变关系如图4-1所示，屈服准则采用莫尔-库仑准则。当土体应力的小于屈服应力的时，不会像刚体一样变形，但一旦到达，即使不变，土体也会像理想塑性体一样发生无限的塑性变形——塑性流动，土体部分或全部进入塑性状态，边坡就会失稳。此时作用在土坡上的荷载等于极限荷载。图4-1刚塑性体的应力应变关系模糊极值理论[64]自1981年以来，石军根据贾政提出的最大最小原理，建立并发展了一种新的稳定性分析理论——模糊极值理论。他首先研究了边坡稳定性分析方法的力学基础，分析了塑性力学极限理论中所包含的模糊性，提出了滑动机制的概念，证明了给定滑动机制的耗散函数定理，并将该定理模糊化，得到了给定滑动面安全系数的模糊极大值定理。然后将合理性要求模糊化为滑动力状态的三个模糊状态条件，构成模糊函数和模糊约束条件。提出了安全系数模糊解集和最小模糊解集的概念，建立了边坡稳定性分析中“最大、中、最小”问题的模糊极值理论。有限元法有限元法是近年来边坡稳定性分析的新趋势。有限元法的优点是可以计算边坡的应力场和位移场的分布，了解边坡的渐进破坏机制，跟踪边坡塑性区的发展。将有限元法应用于渗流计算的数值模拟时，必须首先建立数学模型。描述渗流运动的数学方程，加上初始条件和边界条件，构成了研究特定大坝的数学模型。在数值模拟中，解只是由研究的离散点(大坝渗流场)给出的具有一定精度的未知量的近似值。4.3渗流计算原则数值计算中使用的数学模型按维度不同可分为平面二维渗流和空间三维渗流两种。由于堤坝的长度远大于宽度，所以通常采用平面二维渗流模型。渗流的基本微分方程及其边界条件[65]如图4-2所示；其数学表达式如下:（4-1）初始条件：边界条件：水头边界流量边界式（4-1）为考虑土体压缩的非稳定渗流运动方程，当不考虑土体的压缩性时，，则变为：（4-2）若渗透系数，则式（4-2）即变为典型的拉普拉斯方程。图4-2计算江堤渗流的边界条件4.4渗透稳定判别标准土体抵抗渗透破坏的能力称为抗渗强度，通常以濒临渗透破坏时的水力梯度表示，即临界水力梯度。以土层逸出点渗透比降小于或等于土体允许比降作为安全标准.若,则认为土体产生的渗透变形将失去渗透稳定性,表达式为其中为逸出点的平均渗透比降,分别按有限元方法和文献[66][67]中的公式计算.在有限元计算中,须对渗流场逸出点所在单元的势函数梯度的模求单元平均,即为平均渗透比降。1.流土型土的临界水力梯度[68]当竖向渗流力等于土体的有效重量时，土体就处于流土的临界状态，则流土的临界水力梯度，见式（4-3）（4-3）式中：——土粒比重；——孔隙率从式（4-3）可知，流土临界水力梯度决定于土的物理性质。流土一般发生在渗流的逸出处，因此只需求出渗流逸出处的水力梯度，就可判别流土的可能性：当，则土处于稳定状态；当，则土处于临界状态；当，则土处于流土状态。渗流逸出处的水力梯度实际上是不能求出的，通常把渗流逸出处的流网网格的平均水力梯度作为逸出梯度的。若渗流逸出处网格的水头损失为，网格在流线方向的平均长度为，则逸出梯度为：（4-4）2.管涌型土的临界水力梯度[68]中国水利水电科学研究院根据渗流场中单个土粒受到渗流力、浮力与自重作用时的极限平衡条件，并结合试验资料分析的结果，提出了管涌土临界水力梯度的计算公式为：（4-5）式中：、——小于该粒径的土粒含量分别为5%和20%，其余符号同式（4-3）。综合上述结果，无粘性土的临界水力梯度与其容许值见表4-1表4-1各种土壤的临界水力梯度变化范围及其允许值水力梯度大土集团流动土壤过渡土壤管涌土壤铜>5铜≤5渐变连续性渐变不连续性临界水力陡度0.8~1.01.0~1.50.4~0.80.2~0.40.1~0.3容许水力坡度0.4~0.50.5~0.80.25~0.400.15~0.250.1~0.24.5渗流计算的数值模拟渗流计算方法路堤渗流属于自由面渗流。迄今为止，路堤渗流的计算方法主要有两种，即理论分析法和试验分析法[69-71]。理论分析方法可分为流体力学和水力学。水动力学更严谨，是根据水动力学的基本原理和渗流的边界条件直接求解渗流问题的方法。计算结果更加精确，可以得到渗流场中任意一点的渗流要素(如渗流水头、渗流压力、渗流梯度、渗流速度和通过任意截面的渗流流量)。但这种方法计算复杂，只能用于几种简单的边界条件。水力学方法是基于渗流条件的一些简化假设的解决方案。该方法计算简单，可用于计算各种实际渗流问题。但这种方法只能计算渗流场中某一断面的平均渗流要素，而不能计算渗流场中任意一点的渗流要素。虽然粗糙，但简单实用，精度能满足工程要求，所以应用广泛。根据计算技术的不同，理论分析可以细分为解析法、数值法和图解法。其中，解析法包括直接解法、复变函数法、组合法和水力学法；数值方法分为有限元法和差分法[72][73]；图解法是绘制渗流场中的流网和计算渗流要素的方法[74]。分析方法是指在实验室中用模型模拟实际的水工结构或进行现场观测来研究其渗流要素。主要的模型试验方法有砂槽(土槽)模型法、粘性流模型法(狭缝槽法)、水力网络模型法和水力积分器法、电模拟试验法和电阻网络模型法等。其中，广泛使用的测试方法是模拟测试方法。本文采用美国地质调查局GMS软件中的二维有限元模型SEEP2D进行渗流模拟。通过收集河堤钻孔剖面资料，综合考虑各剖面的工程地质和水文地质性质，选取水位落差最大的断面作为模拟对象，建立了该剖面在恒定流条件下的二维水流模拟模型。4.5.2剖面特征的计算长江大堤主要位于第四纪松散堆积平原。堤基表层相对于隔水层的厚度一般为6.6~8.3m，下部通常为深厚的砂石层。河水和地下水通过沙砾层紧密相连。地基土层的地质结构以二元结构为主，其次为单一结构和多重结构。其中，二元结构地基上部为透水性较弱的粘土、粉质粘土、亚粘土覆盖层，局部夹泥透镜体；下部为透水砂和砂卵石层，透水层厚度可达百米以上，因此该堤基防渗稳定性较差。边界条件的处理河水与其下垫面的交线作为常水头边界处理，下边界作为零流量边界处理。在洪水期，由于地下水位可达地表，岸边的上边界也作为恒定水头边界处理。在模拟过程中，根据长江水位变化调整上下游水头，模拟了不同水头(3、3.5、4、5.5m)下堤下地下水渗流场特征。4.5.4模型划分选取的模型剖面网格为2145个单元，中间粘土层加密。横截面如图4-3所示。图4-3河堤网格剖面图渗流模拟结果根据沉井施工降水情况，选取河堤水位落差最大的三个断面进行渗流分析，编号分别为I-I、II-II、III-III断面。根据模型预测结果，施工期三段河堤上下游水位见表4-2。GMS中SEEP2D模块计算的河岸剖面自由面和等势线见图4-4~图4-8。该图显示了截面II-II的自由面和等势线。表4-2施工期河堤上下游水位(单位:米)第一至第一节第ⅱ-ⅱ节第ⅲ-ⅲ节上游水位下游水位上游水位下游水位上游水位下游水位3.5-0.0583.5-0.0253.5-0.5694.0-0.9264.0-0.9594.0-1.7784.5-1.3334.5-1.2984.5-2.2995.5-1.7485.5-1.8165.5-2.875图4-4沉井不排水河堤的自由面和等势线图4-5沉井下降8m时河堤的自由面和等势线。图4-6沉井下降13m时河堤的自由面和等势线。图4-7沉井下降18m时河堤的自由面和等势线。图4-8沉井降水23m时河堤的自由面和等势线。从等水头线的分布可以看出，地下水流向沉井，流线方向偏向下，水头主要损失在长江固定水头和堤脚之间。从水流速度的分布来看，地面粘土层的流速较小，而下部含水层的流速较大。随着水头差的增大，边界流速逐渐增大。在堤坝渗流问题中，土颗粒级配状况是渗流破坏形成的一个基本条件。土的级配的好坏可由土粒均匀程度和粒径分布曲线的形状来决定，而土粒的均匀程度和粒径分布曲线的形状又可以用不均匀系数（）和曲率系数（）来衡量。由我国《土的分类标准》（GBJ145-90）规定：对于纯净的砾、砂，当大于或等于5，而且等于1～3时，土的级配是良好的，不能同时满足上述条件的，土的级配是不良的。其中，为粒径分布曲线上小于某粒径的土粒含量分别为10%，30%和60%时所对应的粒径。江堤与两层地基土的颗粒分析曲线如图4-9～图4-11。其不均匀系数和曲率系数的计算值如表4-3。从图中可以看出，曲线的坡度较陡，表示土的粒径分布围窄，靠近，土的不均匀系数小，土粒较均匀。根据表4-3中的数据，江堤土和地基土的不均匀系数均小于5，第一层地基土的曲率系数小于1，因此土的级配是不良的，级配不良的土容易发生管涌破坏。因此，在施工降水阶段要对江堤进行监测，并对可能产生的管涌破坏采取预防措施。表4-3和的计算值试样序号江堤土第一层地基土第二层地基土153.22.410.982.080.8220.782.240.893.531.5732.330.861.480.892.330.9644.42.212.410.982.330.9253.251.232.410.986.02.67620.782.290.892.50.961图4-9河堤土颗粒分析试验曲线图4-10第一层地基土颗粒分析试验曲线图4-11第二层地基土颗粒分析试验曲线渗透破坏的另一个原因是水力梯度引起的管涌破坏。当降雨过程中河堤的水力梯度大于土体的临界水力梯度时，土体会发生渗流破坏，威胁河堤的安全。根据降水深度，河堤中心和坡脚产生的水力梯度计算结果见表4-4。表4-4沉井施工降水期间河堤水力坡度河堤水位(米)第一至第一节第ⅱ-ⅱ节第ⅲ-ⅲ节迪申坡脚迪申坡脚迪申坡脚3.50.0480.0760.0590.0710.0690.0624.00.0750.0970.0870.0770.1000.0994.50.0770.0980.0920.0970.1190.1345.50.2450.130.2480.1360.1430.157从上表可以看出，河堤中心与坡脚之间的水力梯度随着上下游水位差的增大而增大。在最大水位差作用下，水力比降较大，特别是在河堤中心，接近临界水力比降。表4-4中计算的水力梯度都在表4-1中列出的允许水力梯度范围内。因此，从水力坡度的判据来看，在沉井施工降水过程中，堤防不会受到渗流的破坏。4.6边坡稳定性分析与计算计算方法大堤的地质结构主要是二元结构。在稳定性计算中，考虑了一些特殊的外力，如渗透力和地震惯性力。因此，选用瑞典圆弧滑动法来验算其抗滑稳定性。瑞典圆弧滑动法(简称瑞典法)假设滑动面为圆柱面(截面上为圆弧)，同时也假设不考虑土条两侧的力。安全系数定义为滑动面上每个土条提供的抗滑力矩之和与滑动面上滑动土产生的外荷载和滑动力矩之和的比值。应力分析如图4-12和4-13所示。考虑到渗透力，瑞典方法的计算公式为:式中：——土的有效抗剪强度指标；——水的容重。图4-12瑞典法计算图图4-13渗流对土坡稳定的影响4.6.2计算断面与条件抗滑稳定的计算断面选用渗流计算的断面一致，如图4-2所示。江堤外侧的水位根据预测的长江水位变化，采用3.0m、3.5m、4.0m和5.5m四种水位。4.6.3计算结果根据《堤防设计规[75]》，江堤的抗滑稳定安全系数不应小于表4-5的规定。本文采用瑞典法计算的抗滑稳定安全系数见表4-6。Ⅱ-Ⅱ剖面不同降水时段的稳定计算如图4-14～图4-17。表4-5江堤的抗滑稳定安全系数堤防工程的级别12345安全系数正常运用条件1.301.251.201.151.10非常运用条件1.201.151.101.051.05表4-6江堤抗滑稳定安全系数水位（m）沉井降深（m）Ⅰ-Ⅰ剖面Ⅱ-Ⅱ剖面Ⅲ-Ⅲ剖面3.581.5901.5881.6024.0131.6131.6311.6514.5181.6271.6261.6515.5231.6071.6081.575图4-14沉井降水8m图4-15沉井降水13m沉井降水18图4-17沉井降水23m注:图中虚线表示外堤水位线和堤防浸润线，地面标高为2.5m从图表中可以看出，在施工降水过程中，抗滑稳定安全系数大于河堤正常运行工况下的抗滑稳定容许安全系数。因此，堤防的抗滑稳定是安全的。为了判断河堤抗滑安全系数随上下游水位差的变化趋势，根据表4-2和表4-6绘制的水位差与抗滑安全系数的关系曲线如图4-18所示。图4-18水位差与抗滑稳定安全系数关系曲线从图4-18可以看出，上下游水位差较小时，抗滑稳定系数随着上下游水位差的增大而增大；当上下游水位差继续增大时，抗滑稳定系数随着上下游水位差的增大而减小。这一变化规律表明，在水位下降过程中，边坡的稳定系数有一个最大值，造成这一结果的原因是边坡中的浸润线也发生了变化。随着水位的降低，稳定系数按照由大到小的规律变化，与参考文献[76]提出的规律一致，说明计算是合理的。水位下降初期，稳定系数随着上下游水位差的增大而增大，因为浸润线的降低增加了土的有效应力和土的自重，增加的向下抗滑力大于水平向外渗流引起的土的滑动力。然后，随着水位差的进一步增大，河堤的水力梯度逐渐增大，渗透力引起的土体下滑力大于土体自身重量增大的向下抗滑力。因此，抗滑稳定系数逐渐减小。4.7本章概述本章首先详细阐述了渗流理论的破坏形式、计算方法和稳定判据，并用GMS软件中的二维有限元模型SEEP2D进行模拟，得到完整的流网和流线，分析了河堤在渗透力作用下的渗透稳定性。然后，根据渗流计算的自由面位置，采用瑞典圆弧滑动法计算河堤抗滑稳定安全系数，得到抗滑稳定安全系数随上下游水位差增大先增大后逐渐减小的变化规律。6主要结论和展望6.1本文的主要工作和成果本文结合在建长江大桥超大型沉井施工降水方案研究，利用地下水模拟软件GMS建立水源地数值模型，评价施工期地下水涌水量，预测施工期地下水动态。通过研究，得出以下结论:1.本文将研究区域概化为一个非均质、各向异性、三维结构和不稳定的地下水流系统，首次将GMS应用于特大沉井施工降水分析。本研究采用概念模型和网格模型进行建模，简单方便，提高了建模的灵活性。GMS的高度可视化特性可以更好地再现地下水的动态变化。2.根据地下水数值模型的识别和验证标准，对水平和垂直渗透系数进行反复调整和识别，将抽水量、静态地下水位和水位降落深度的实测和统计数据代入模拟模型，得到理想的拟合流场。模型验证后，实际观测水位与计算水位之差小于0.5m3.模型采用施工降水方案进行预测，模型预测的涌水量与一般理论公式计算的涌水量基本一致。这进一步验证了模型的可靠性和准确性，为施工的顺利进行提供了可靠的依据。从等水头图可以看出，在施工阶段，随着水位的降低，网格排水的深度和周长逐渐增大，形成的降落漏斗水力梯度由外向内由小变大，抽水井区附近的水力梯度最大，符合渗流规律。从降水过程线的变化趋势可以看出，靠近长江的沉井两侧水位下降缓慢，周长较小。这是由于扬子江在抽水过程中不断补给，而其他两岸水量不足，所以水位下降快，周长大。4.从渗流基本微分方程出发，利用二维有限元模块SEEP2D建立模型，求解研究区地下水的自由面、流网和流线，得到完整的流网和流线。同时，分析了水位下降时边坡失稳机理和抗滑安全系数的变化规律。结果表明，在水位下降过程中，当上下游水位差较小时，堤防抗滑稳定系数随着上下游水位差的增大而增大，随着上下游水位差的增大而减小。5.对于地面沉降的计算，本文考虑了抽水后含水层的垂直主固结。用BCF程序对平面问题的Biot固结方程进行有限元求解，得到了沉降随时间的变化。从计算结果可以看出，最大沉降发生在河堤坡脚附近。根据《堤防工程设计规范》(GB50286-98)，本研究区的沉降在理论控制范围内。6.2论文的创新点1.利用真实的地层数据进行模拟，可以对模型进行反复修正，使其与现场观测值一致。首次将GMS软件应用于特大沉井施工降水分析。2.本课题采用三维可视化动态模拟，可以输出动态水流路径，获得完整的流网和等势线，准确分析河堤的渗透稳定性。3.现行规范尚未提出大坝的允许沉降圈。该模拟可以预测大坝最大沉降段的位置和最大沉降量，为控制河堤沉降变形提供依据。4.长江公路大桥是国家重点工程。本课题将以模拟结果指导现场施工方案的实施，进一步验证模型计算的准确性。整个沉井施工正按计划顺利进行。6.3有待进一步研究的问题和展望1.由于研究