地下水数值模拟-第三讲-地下水数值模拟原理及建模方法和步骤课件.ppt

地下水数值模拟原理及建模方法和步骤,内容提要,绪论一、地下水流有限差分法原理二、地下水数值模型建模步骤三、建模所需要的基本资料 四、地下水系统溶质迁移有限差分法原理五、地下水系统热对流-扩散有限差分法原理,绪论地下水数值方法在水文地质学中的位置,已学课程水文地质学基础地下水动力学水文地球化学（环境水文地质学）同位素地下水溶质运移：示踪剂或污染物在含水层中的迁移机理及数学模型和求解方法地下水数值模拟,绪论地下水数值方法在水文地质学中的位置,地下水动力学主要内容连续性原理、达西定律、水均衡原理地下水流基本方程几类特殊水文地质问题数学模型及解析解地下水向沟渠河中的流动园岛模型泰斯模型有越流的不稳定井流（Hantush and Jacob）无越流的潜水含水层不稳定井流（Neuman）抽水试验及反求参数数值方法可以作为地下水动力学课程内容之一、或补充、或延续,绪论数值解与解析解,地下水动力学中所得到的解是解析解解析解的特点可以利用解析解公式计算出研究区任意空间点和时间点处的水头值对实际模型的概化适用于一些特殊简单问题公式复杂，需要借助计算机求解因此，其应用受到很大限制,绪论数值解与解析解,数值方法是地下水动力学的完善和补充或延续数值解的特点：只是求出研究区某些空间点和某些时间点处的水头值适用于所有的问题具备水文地质基础和线性代数知识已有数值模拟专门软件（或自己编程）需要有高性能计算机对实际问题的刻画比较精确因此，其应用非常广泛,数值方法很多，但是最简单实用的是有限差分法：,绪论数值方法,有限差分法有限单元法积分有限差分法半解析半数值法边界元法有限体积法,只讲有限差分法,一、有限差分法的基本原理,有两种方法建立差分方程方法一以地下水流基本微分方程及其定解条件为基础，在渗流区剖分基础上，用差商代替微商，将地下水流微分方程的求解转化为差分方程（代数方程）求解。适用于二维矩形网格剖分、三维长方体网格剖分。方法二在渗流区剖分的基础上，直接由达西定律和水均衡原理，建立各个均衡区的水均衡方程，从而得到差分方程。适用于矩形网格、三角形网格。,矩形网格,多边形网格,1、网格划分的基本类型,（1）先划格线，格点位于网格中心,均衡网格,节点网格,（2）先规定格点位置，再垂直平分两相邻结点的连线作格线，形成的网格即为水均衡区,方法一：差商代替微商,MODFLOW网格系统,方法一：差商代替微商,导数的有限差商近似,导数的定义,当,非常小的时候，有,上式右端项即为f(x)在x0处的差商。,这样定义的差商很容易理解，但不知道用差商代替微商所产生的误差。下面利用泰勒公式导出差商及其误差。,方法一：差商代替微商,(2)有限差分方程建立,已知泰勒公式,由A得：,A,B,由B 得：,称 为f(x)在x0处的一阶前向差商，为截断误差。,称 为f(x)在x0处的一阶后向差商，为截断误差。,方法一,由A-B可以得：,由A+B可以得：,A,B,称 为f(x)在x0处的一阶中心差商，为截断误差。,称 为f(x)在x0处的二阶中心差商，为截断误差。,方法一,对于偏导数（偏微商），类似可以得到相应的差商：,方法一,(2)有限差分方程建立(续),一维控制方程差分格式,显式差分格式,隐式差分格式,方法一,控制方程,网格剖分nx个,二维控制方程差分格式,显式差分格式,隐式差分格式,方法一,控制方程,网格剖分nx个,达西定律：水均衡原理：对某一研究对象，流入 流出体系内质量（或水量）变化量研究对象可以是大区域的，也可以是微分单元体大区域的水均衡计算经常用于区域的水资源评价,单位时间通过单位面积过水断面的水量与断面处的水力梯度成正比。,方法二：直接由达西定律和水均衡原理建立差分方程,表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）,+,H 可取tn+1或 tn时刻的值,节点（i,j）的均衡区,Aij:节点(i,j)均衡区的面积,不等距矩形网格有限差分法：以二维承压水流为例,表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）,H取tn时刻的值，有：,+,节点（i,j）均衡区,Aij:节点(i,j)均衡区的面积,不等距矩形网格有限差分法：以二维承压水流为例（续1）,则得到显式格式：,不等距矩形网格有限差分法：以二维承压水流为例（续2）,+,如果 H 取tn+1时刻的值,节点（i,j）的均衡区,不等距矩形网格有限差分法：以二维承压水流为例（续3）,记：,得到隐式格式：,不等距矩形网格有限差分法：以二维承压水流为例（续4）,表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）,当网格是等距时，即,+,H 可取tn+1或 tn时刻的值,节点（i,j）的均衡区,此时，有水均衡方程：,不等距矩形网格有限差分法：以二维承压水流为例（续5）,节点（i,j）的均衡区,两边除以，得到,表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）,H 取tn+1时刻的值，得到隐式格式,不等距矩形网格有限差分法：以二维承压水流为例（续6）,表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）,两边除以,+,H 可取tn+1时刻的值，得到隐式格式,不等距矩形网格有限差分法：以二维承压水流为例（续7）,表示：当水头下降一个单位时，从单位底面积饱水空隙介质中释放的水量（体积）,如果（i,j)是左边界节点，有,+,H 可取tn+1或 tn时刻的值,节点（i,j）的均衡区,qij:节点(i,j)处的边界单宽流量。,不等距矩形网格有限差分法：以二维承压水流为例（续8）,三维有限差分格式（显示）,三维有限差分格式（隐式）,首先将渗流区划分成若干个辅助小三角形（图36中虚线所示）。划分时应注意：三角形的任一内角不得大于90，三条边的长度尽可能接近；三角形的顶点不能落在另外某个三角形的边上；相邻三角形变化不太大；应考虑水文地质条件，尽量使抽水井、观测孔（特别是用于拟合的观测孔）位于三角形的顶点处，河流、断层等位于三角形边上。,1渗流区的剖分,任意多边形网格系统,三角形的顶点称为格点，分为内格点和边界格点（包括第一类边界格点和第二类边界格点）。,2、任意多边形网格有限差分法,一、有限差分法的基本原理,2多边形均衡区网格的形成方法是：对某个格点i，分别作出格点i与其相邻格点的连线的垂直平分线，这些平分线所围成的多边形区域即为格点 i 的多边形均衡网格子区,对每个格点都这样做，则形成多边形均衡网格系统（图36中实线所示）,典型多边形子区,2、任意多边形网格有限差分法,典型多边形子区,i=1,2,N,隐式格式,水均衡方程为：,2、任意多边形网格有限差分法,典型多边形子区,i=1,2,N,显式格式,2、任意多边形网格有限差分法,（1）区域剖分 首先，在平面上将研究区剖分成若干个辅助小三角形，然后以这些三角形为底面，将区域剖分成一些垂直的柱形体，再根据地层岩性及构造特征，用一些平面将柱形体剖分成若干短柱体（图3.1）。关于平面三角形的剖分依据下述原则：1）三角形的任一内角不得大于90，三条边的长度尽可能接近；2）三角形的顶点不能落在某个三角形的边上；3）相邻三角形变化不太大。,3、多棱柱体三维流有限差分法,（2）均衡区网格的形成 三角形的顶点称为格点，分为内格点和边界格点（包括第一类边界格点和第二类边界格点）。各结点的均衡区形成方法是：平面上，以格点为中心，与之相连的边的垂直平分线段组成一个多边形；垂向上，以该层与上下两层的两个中间层的平面分别为上下底面，这样得到一个以平面上为多边形的立体柱体，就是水均衡子区。,3、多棱柱体三维流有限差分法,考虑以节点（i,m）为中心的多面柱体均衡子区，下面根据达西定律和水均衡原理建立其差分方程。如图所示，单位时间内通过单元中两断面pb、bq 断面流入到均衡区的水量为:,(i,m),3、多棱柱体三维流有限差分法,式中：,所对应的柱面流入均衡区内的水量；,为第m层流段ij和ik单元的导水系数，定义为,从第e号三角形通过两线段,流段的长度。,类似对格点 i 周围所有三角形作上述计算，并求和得到从侧面流入到均衡区内的总水量,为,另外，容易求得通过均衡区顶底面流入到均衡区内的水量为,式中：,第i个均衡区的平面投影区域面积；,第i个平面格点对应第m层和m-1层之间z方向的平均渗透系数；,第i个平面格点对应第m层和m-1层之间z方向的平均渗透系数；,第i个平面格点对应第 m 层和 m-1层之间z方向的平均渗透系数；,若记 为第i个平面格点对应第m层的均衡子区的源汇项，则由水均衡原理可获得三维地下水差分方程：,式中：,N研究区平面网格结点总数，,Mi第i个平面结点对应的垂向分层数。,比储水系数,对于非稳定流动问题，由初始时刻的水头分布开始，利用上时段的水头值递推计算以后各个时段的水头值，因此公式中为已知，为待求水头，将方程整理得：,-,4、潜水（无压）水流有限差分法,潜水含水层平面二维流,潜水含水层剖面二维流,对于剖面二维流，完全可以作为三维流的一种特例处理，因此这里只讨论平面二维流情形。,4.1显-隐式法,潜水含水层与承压含水层不同，随着潜水面的升降，会引起导水系数的变化，所以潜水含水层的导水系数是时间的函数。所谓凝固系数法就是将某时段的潜水含水层的导水系数利用该时段初的水头值计算，即令,从而，潜水流方程可写成,4、潜水（无压）水流有限差分法,3.7.2 全隐式法,潜水含水层与承压含水层不同，随着潜水面的升降，会引起导水系数的变化，所以潜水含水层的导水系数是时间的函数。所谓双重迭代法就是利用迭代法求解非线性潜水流微分方程。在第m迭代步，潜水含水层的导水系数利用该迭代步的水头值计算，即令,求解方程,得到H（m），直到收敛为止。,4、潜水（无压）水流有限差分法,3-8-1 越流、入渗和抽水井等问题的处理如果考虑垂直渗流项（即源汇项），则二维承压流动微分方程可写成建立差分方程时，在结点处应加上这一项，它可具体表示为式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相邻弱透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含水层水头；为主含水层水头；为入渗强度或蒸发强度；为抽水或注水流量。,5、特殊问题的处理：井孔水头校正,井水位校正对比图,图3-5（a）大步长网格 图3-5（b）加密网格,有限差分法计算井水位的校正,Q0 为注水；Q0 为抽水,5、特殊问题的处理：不规则边界问题,当研究区的几何形状属于简单形式（如矩形渗流区）时，差分网格的划分往往将结点设在边界上。然而，对于实际问题来说，边界通常不是那么规则，边界的某些部分，甚至大部分不能与结点重合，我们称这种边界是不规则的。关于不规则边界问题，直接取最靠近边界的曲折格线为近似边界,越流、入渗和抽水井等问题的处理如果考虑垂直渗流项（即源汇项），则二维承压流动微分方程可写成建立差分方程时，在结点处应加上这一项，它可具体表示为式中：为垂向渗流强度（包括越流、入渗或蒸发、井流等）；为相邻弱透水层垂向渗透系数；为相邻弱透水层厚度；为相邻含水层水头；为主含水层水头；为入渗强度或蒸发强度；为抽水或注水流量。,源汇项的处理:井孔水头校正,6、差分方程求解,一维显式差分格式,网格个数为ni,直接求解,6、差分方程求解,一维隐式差分格式,网格个数为ni,迭代求解,方程组,PCGSIPSORWHSSAMGGMG,MODFLOW,7、差分方程的收敛性和稳定性,截断误差：用差商代替微商时，地下水流动方程产生的误差为截断误差。收敛性：当空间步长和时间步长趋于0时，有限差分方程的精确解趋于地下水流动问题微分方程定解问题的精确解。则称该差分格式是收敛的。稳定性：如果在求解差分方程过程中，某时间步引入某个误差，而在以后的各时段计算中，该误差不再扩大，则称该差分格式是稳定的。,一维显示格式的收敛条件和稳定条件是：,8、算例：显式有限差格式,设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。,应用实例：河间地块承压水流模型,步骤：,（1）基础资料的分析（2）概念模型（3）数学模型（4）数值方法及计算机程序（5）参数（6）结果分析,8、算例：显式有限差格式（续1）,建立数学模型,（1）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。（2）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内，从而量化各变量的取值范围。本例，取x-轴原点位于左端河，右侧为正向，设两河流间距为L.（3）数学模型,8、算例：显式有限差格式（续2）,将（0L）分成 N 等份，,1)网格剖分：,取时间步长，记（n=0、1、2、3、4）,记，(i=0,1,2,3，4N),2）建立差分方程：在网格系统中任意取一点,设,是问题的解，则在,处有,记为（i，n）,8、算例：显式有限差格式（续3）,用差商代替微商：,将上述两式舍去余项，代入方程并记,为,显然该式具有截断误差,得到,8、算例：显式有限差格式（续4）,引入无量纲变量：,将该式子代入得到：,（i=1，2，3，.N-1），（n=1，2，3，.）,8、算例：显式有限差格式（续5）,3）显示差分方程的求解,计算各结点初始时刻水头值利用差分方程计算各结点t1时刻水头值利用边界条件计算边界结点水头值重复2、3步，直到计算出拟计算的各个时刻的水头值,8、算例：显式有限差格式（续6）,在上述模型中，设L=1000米,取空间步长为200米，时间步长为0.25天，分别计算各节点各时刻的水头值。,8、算例：显式有限差格式（续7）,8、算例：显式有限差格式（续8）,如果 t=1,则,8、算例：显式有限差格式（续9）,在上述模型中，设L=1000米,取空间步长为200米，时间步长为0.25天，用隐式差分格式计算各节点个时刻的水头值。,8、算例：隐式有限差格式,在这个例子中，,解：隐式格式一般方程为,于是有,根据初始条件得,根据边界条件得,由初始条件和边界条件,由此解得t1时刻的水头值为,在上述方程中取 n=0,可以得到计算t1时刻水头值的方程,所以上述方程变成,同理，可计算t2时刻的水头值,给定一组试探解,分别代入上述三个方程，计算出,简单迭代法,对于方程组,如果 是原方程组的解，则 与 非常接近。,8、算例：显式有限差格式（续2）,通常给定一个很小的正数，若,在计算出第K步的迭代值后，按下式计算第K+1步的迭代值，,则称 为原方程组的近似解。,一般情形，不可能一次就找到方程组的近似解，通常要计算很多步才可能得到满足精度要求的近似解。,若 则取 为原方程组的近似解。,8、算例：显式有限差格式（续3）,简单迭代法的收敛条件：如果方程组的系数矩阵是对角占优矩阵，则相应的迭代格式是收敛的。由收敛迭代格式计算出的向量序列 一定收敛到某个确定的向量，该向量就是原方程组的精确解。,8、算例：显式有限差格式（续1）,迭代法-改进迭代法,在计算出第K步的迭代值后，按下式计算第K+1步的迭代值，,一般情形，对于一个收敛的迭代格式来说，第k+1步的迭代值会比第K步的迭代值更接近精确解，因此，根据这一认识，可以对简单迭代格式进行改进。,若 则取 为原方程组的近似解。,迭代法-改进迭代法（续）,通常给定一个很小的正数，若,在计算出第K步的迭代值后，按下式计算第K+1步的迭代值，,则称 为原方程组的近似解。,一般情形，不可能一次就找到方程组的近似解，通常要计算很多步才可能得到满足精度要求的近似解。,若 则取 为原方程组的近似解。,内容提要,绪论一、地下水流有限差分法原理二、地下水数值模型建模步骤三、建模所需要的基本资料 四、地下水系统溶质迁移有限差分法原理五、地下水系统热对流-扩散有限差分法原理,二、地下水数值模型建模步骤,模拟步骤,建立概念模型 建立数学模型 数值方法及软件（编程）参数准备以及计算 模型校正与检验 参数敏感性分析 预测模拟,软件,一、概念模型（模型概化）,根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、岩石矿物、水文、气象、工农业利用情况等模拟的区域:含水层类型:潜水（无压）、承压、混合、多层维数:一维、二维、三维水流状态：稳定流/非稳定流、饱和流/非饱和流介质状况:均质和非均质/各向同性和各向异性孔隙/裂隙/双重介质流体的密度差边界条件和初始条件 必要时需进行一系列的室内试验与野外试验,以获取有关参数,如渗透系数、弥散系数、分配系数、反应速率常数等。,（二）数学模型,三维地下水流动问题控制方程,第二类边界条件,第一类边界条件,初始条件,绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。,三、数值方法及软件（或编程）,有限差分法有限单元法积分有限差分法半解析半数值法边界元法,有限差分法:MOFLOW系列GMS 中MODFLOWVisual MODFLOWProcessing MODFLOW有限单元法:FEFLOW积分有限差分法:TOUGH2，TOUGH REACT,三、数值方法及软件（或编程）,四、模型参数准备以及计算,含水层参数：渗透系数，弹性释水系数，孔隙度等源汇项：大气降水入渗系数（分区、数值）蒸发排泄系数地表水体水位、底面高程、底面岩性特征渠系灌溉入渗系数人工开采（点状、面状）边界条件初始条件,计算,资料,结果,输入模块,运行模块,输出模块,Visual MODFLOW,参数不确定性,钻孔太少，地层资料少，钻孔多，含水层结构会发生变化,四、模型校正与检验,（1）钻孔资料（2）抽水试验,参数不确定性,泰斯井流抽水试验假设：（1）（2）（3）（4）（5）（6）,（四）模型校正与检验,（1）钻孔资料（2）抽水试验,（五）模拟：模型校正（参数识别）,将模拟结果与实测结果比较,进行参数调整,使模拟结果在给定的误差范围内与实测结果吻合。调参过程是一个复杂而辛苦的工作,所调整的参数必须符合模拟区的具体情况。人机交互与自动调参相结合。尽管自动调参程序(如PEST),也不能代替人的工作。,（五）模拟：模型检验,模型验证是在模型校正的基础上,进一步调整参数,使模拟结果与第二次实测结果吻合,以进一步提高模型的置信度。,（六）灵敏度分析,校正后的模型受参数值的时空分布、边界条件、水流状态等不确定度的影响。灵敏度分析就是为了确定不确定度对校正模型的影响程度。,（七）预测,用校正的参数值进行预测,预测时需估算未来的水流状态。,后续检查与模型的重建（完善）,后续检查在模拟研究结束数年后进行。收集新的野外数据以确定预测结果是否正确。如果模拟结果精确,则该模型对该模拟区来说是有效的。由于场址的唯一性,故模型只对该模拟区有效。后续检查应在预测结束足够长的时间后进行,以便有足够的时间发生明显的变化。,模型的再设计,一般来说,后续检查会发现系统性能的变化,从而导致概念模型和模型参数的修改。一般来说,所有模拟研究都应该进行到第五步,即校正灵敏度分析。,地裂缝,自然地理及水文地质条件,边界敏感性分析,数学模型,网格化,概念模型,介质类型、结构特征,地下水补、径、排特征,数值模型,数值模型的建立,边界条件,初始条件,源汇项,降雨入渗,地表水体入渗,灌溉、渠系,蒸发排泄,人工开采,介质参数,观测孔动态拟合,流场拟合,水均衡对比分析,观测孔动态检验,流场检验,水均衡对比检验,拟合调参,识别结果不符合要求修正概念模型,小结：地下水流动模型构建过程,识别结果不符合要求修正概念模型,模型应用,水资源量评价,预测开采方案,研究环境地质问题,地面沉降,地下水污染研究,大型工程对环境的影响,模型识别,模型检验,内容提要,绪论一、地下水流有限差分法原理二、地下水数值模型建模步骤三、建模所需要的基本资料 四、地下水系统溶质迁移有限差分法原理五、地下水系统热对流-扩散有限差分法原理,三、地下水流数值模拟资料需求,空间展布相关资料地表高程、分层数据边界位置含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等源汇项：大气降水入渗系数（分区、数值）蒸发排泄系数地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等）渠系灌溉入渗系数人工开采（点状、面状）：开采井位置、井结构、开采量动态边界条件：边界类型、水头或流量初始条件：统测水位水位动态观测资料：观测孔位置、结构、水位时间变化,内容提要,绪论一、地下水流有限差分法原理二、地下水数值模型建模步骤三、建模所需要的基本资料 四、地下水系统溶质迁移有限差分法原理五、地下水系统热对流-扩散有限差分法原理,在多组分组成的流体体系中任取一点P(x,y,z),以P为中心取一微小的质量平衡体，其侧面分别平行与3个坐标面，边长分别为x、y、z。,质量守恒原理:在时间t内，组分在这个单元体中的净流出（或流出）量（暂不考虑起内部有质量产生和消失），应等于这个单元中组分的质量变化,1、污染物迁移有限差分方程建立,X+Y+Z+I=m,Fick定律和质量守恒原理,1、污染物迁移有限差分方程建立(续),这个单元中组分的质量变化量,有质量守恒原理X+Y+Z+I=m，可得：,隐式差分格式,n+1,显式差分格式,n,内容提要,绪论一、地下水流有限差分法原理二、地下水数值模型建模步骤三、建模所需要的基本资料 四、地下水系统溶质迁移有限差分法原理五、地下水系统热对流-扩散有限差分法原理,在多组分组成的流体体系中任取一点P(x,y,z),以P为中心取一微小的能量平衡体，其侧面分别平行与3个坐标面，边长分别为x、y、z，,能量守恒原理:在时间t内，在这个单元体中的净流出（或流出）量+内部有能量产生和消失，应等于这个单元中能量变化,1、热迁移有限差分方程建立,X+Y+Z+I=m,Fourier定律和能量守恒原理,断面流入的能量由两部分构成：对流作用与弥散作用，其中对流作用采用Fourier定律来刻画，x方向变化量：,弥散：,2.3热迁移有限差分方程建立(续),对流：,X,X轴方向能量变化量,+,+,弥散：,2.3热迁移有限差分方程建立(续),对流：,Y轴方向能量变化量,+,+,弥散：,对流：,Z轴方向能量变化量,+,+,Y,Z,1、热迁移有限差分方程建立(续),这个单元中能量变化量,有质量守恒原理X+Y+Z+I=m，可得：,隐式差分格式,n+1,显式差分格式,n,小结,一、地下水流有限差分法原理二、地下水数值模型建模步骤三、建模所需要的基本资料 四、地下水系统溶质迁移有限差分法原理五、地下水系统热对流-扩散有限差分法原理,