第二讲 地下水建模方法和步骤精选文档.ppt

第二讲 地下水建模方法和步骤本讲稿第一页，共四十九页地下水建模方法和步骤地下水建模方法和步骤1.求解地下水运动方程的数值方法求解地下水运动方程的数值方法2.地下水数值模型建模步骤地下水数值模型建模步骤3.建模所需要的基本资料建模所需要的基本资料 本讲稿第二页，共四十九页 绝大部分数学模型是无法用解析法求解的，数值化就绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。是将数学模型转化为可解的数值模型。1.数值方法数值方法1.1有限差分法有限差分法1.2有限单元法有限单元法1.3分有限差分法分有限差分法1.4半解析半数值法半解析半数值法1.5边界元法边界元法本讲稿第三页，共四十九页o(1)有限差分法原理有限差分法原理o(2)两种方法建立有限差分方程两种方法建立有限差分方程o(3)求解有限差分方程求解有限差分方程o(4)收敛性和稳定性概念收敛性和稳定性概念o(5)算例算例1.1有限差分法有限差分法本讲稿第四页，共四十九页(1)有限差分法的基本原理有限差分法的基本原理将连续的问题离散后求解：将连续的问题离散后求解：方法一方法一以地下水流基本微分方程以地下水流基本微分方程及其定解条件为基础，及其定解条件为基础，在渗流区剖在渗流区剖分基础上，用差商代替微商，将地下分基础上，用差商代替微商，将地下水流微分方程的求解转化为差分方程水流微分方程的求解转化为差分方程（代数方程）求解。（代数方程）求解。方法二方法二在渗流区剖分的基础上，直接在渗流区剖分的基础上，直接由达西定律和水均衡原理，建立各个均由达西定律和水均衡原理，建立各个均衡区的水均衡方程，即差分方程。衡区的水均衡方程，即差分方程。矩形网格矩形网格多边形网格多边形网格本讲稿第五页，共四十九页网格划分的基本类型网格划分的基本类型o（1）先划格线，格）先划格线，格点位于网格中心点位于网格中心均均衡衡网网格格节节点点网网格格o（2）先规定格点位置，）先规定格点位置，再垂直平分两相邻结点的再垂直平分两相邻结点的连线作格线，形成的网格连线作格线，形成的网格即为水均衡区即为水均衡区本讲稿第六页，共四十九页MODFLOW网格系统网格系统本讲稿第七页，共四十九页导导数数的的有有限限差差商商近近似似导数的定义导数的定义 当当非常小的时候，有非常小的时候，有 上式右端项即为上式右端项即为f(x)(x)在在x0 0处的差商。处的差商。这样定义的差商很容易理解，但不知道用差商代替微商所产生的误差。下面利用泰勒公式导出差商及其误差。方法一：方法一：差商代替微商差商代替微商(2)有限差分方程建立有限差分方程建立本讲稿第八页，共四十九页已知泰勒公式 由A得：AB 由B 得：称称 为为f(x)在在x0处的处的一阶前向差商，一阶前向差商，为为截断误差截断误差。称称 为为f(x)在在x0处的处的一阶后向差商，一阶后向差商，为为截断误差截断误差。方法一方法一本讲稿第九页，共四十九页 由A-B可以得：由A+B可以得：AB称称 为为f(x)在在x0处的处的一阶中心差商，一阶中心差商，为为截断误差截断误差。称称 为为f(x)在在x0处的处的二阶二阶中心中心差商，差商，为为截断误差截断误差。方法一方法一本讲稿第十页，共四十九页l对于偏导数（偏微商），类似可以得到相应的差商：对于偏导数（偏微商），类似可以得到相应的差商：方法一方法一(2)有限差分方程建立有限差分方程建立(续续)本讲稿第十一页，共四十九页一维控制方程差分格式一维控制方程差分格式显式差分格式显式差分格式隐式差分格式隐式差分格式方法一方法一控制方程控制方程网格剖分网格剖分nx个个本讲稿第十二页，共四十九页 取右图所示得微小六面取右图所示得微小六面体。设与体。设与x,y,z,方向对应得主方向对应得主渗透系数分别为渗透系数分别为Kx,Ky,Kz；建立均衡期；建立均衡期 t时段内，微时段内，微小均衡六面体的水量守恒小均衡六面体的水量守恒方程。方程。方法二：方法二：达西定律和水均衡原理达西定律和水均衡原理(2)有限差分方程建立有限差分方程建立(续续)本讲稿第十三页，共四十九页基于基于达西定律达西定律，x，y，z方向流入方向流入流出分别为流出分别为:l t时段内，侧向流入与源汇项导致六面体水量变化量为：时段内，侧向流入与源汇项导致六面体水量变化量为：ABC(2)有限差分方程建立有限差分方程建立(续续)方法二：方法二：达西定律和水均衡原理达西定律和水均衡原理DA+B+C+D源汇项源汇项本讲稿第十四页，共四十九页六面体内地下水储存量的变化为六面体内地下水储存量的变化为由水均衡原理得三维地下水流动方程的有限差分格式由水均衡原理得三维地下水流动方程的有限差分格式(2)有限差分方程建立有限差分方程建立(续续)方法二：方法二：达西定律和水均衡原理达西定律和水均衡原理本讲稿第十五页，共四十九页有限差分法有限差分法:三维（三维（MODFLOW）差商代替微商差商代替微商本讲稿第十六页，共四十九页(3)差分方程求解差分方程求解l一维显式差分格式网格个数为网格个数为ni直接求解直接求解本讲稿第十七页，共四十九页(3)差分方程求解差分方程求解l一维隐式差分格式一维隐式差分格式网格个数为网格个数为ni迭代求解迭代求解方方程程组组PCGSIPSORWHSSAMGGMGMODFLOW本讲稿第十八页，共四十九页(4)差分方程的收敛性和稳定性差分方程的收敛性和稳定性l截断误差截断误差：用差商代替微商时，地下水流动方程产生的误用差商代替微商时，地下水流动方程产生的误差为截断误差。差为截断误差。l收敛性收敛性：当空间步长和时间步长趋于当空间步长和时间步长趋于0 0时，有限差分方程的时，有限差分方程的精确解趋于地下水流动问题微分方程定解问题的精确解。精确解趋于地下水流动问题微分方程定解问题的精确解。则称该差分格式是收敛的。则称该差分格式是收敛的。l稳定性稳定性：如果在求解差分方程过程中，某时间步引入某如果在求解差分方程过程中，某时间步引入某个误差，而在以后的各时段计算中，该误差不再扩大，个误差，而在以后的各时段计算中，该误差不再扩大，则称该差分格式是稳定的。则称该差分格式是稳定的。一维一维显示格式显示格式的收敛条件和稳定条件是：的收敛条件和稳定条件是：本讲稿第十九页，共四十九页(6)算例：算例：显式有限差格式显式有限差格式 设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。应用实例：河间地块承压水流模型应用实例：河间地块承压水流模型本讲稿第二十页，共四十九页步骤：步骤：（1 1）基础资料的分析）基础资料的分析（2 2）概念模型）概念模型（3 3）数学模型）数学模型（4 4）数值方法及计算机程序）数值方法及计算机程序（5 5）参数）参数（6 6）结果分析）结果分析 本讲稿第二十一页，共四十九页建立数学模型建立数学模型（1 1）模型概化）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。由所述水文地质条件，可以概化为一维承压水流问题。（2 2）建立坐标系）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内，从而量（如图），将地下水流动系统空间结构放在坐标系内，从而量化各变量的取值范围。本例，取化各变量的取值范围。本例，取x-x-轴原点位于左端河，右侧为正向，设轴原点位于左端河，右侧为正向，设两河流间距为两河流间距为L.L.（3 3）数学模型）数学模型本讲稿第二十二页，共四十九页差分方程及其解法差分方程及其解法显式格式显式格式将（0L）分成 N 等份，1)1)网格剖分：网格剖分：取时间步长取时间步长 ，记，记 （n n=0=0、1 1、2 2、3 3、44）记 ，(i=0,1,2,3，4N)2 2）建立差分方程）建立差分方程：在网格系统中任意取一点：在网格系统中任意取一点设设是问题的解，则在是问题的解，则在处有处有记为（记为（i i，n n）本讲稿第二十三页，共四十九页用差商代替微商用差商代替微商：将上述两式舍去余项，代入方程并记将上述两式舍去余项，代入方程并记为为显然该式具有截断误差得到得到显式格式（续显式格式（续1）本讲稿第二十四页，共四十九页引入无量纲变量引入无量纲变量：将该式子代入得到：将该式子代入得到：（i=1，2，3，.N-1），（n=1，2，3，.）显式格式（续显式格式（续2）本讲稿第二十五页，共四十九页3）显示差分方程的求解）显示差分方程的求解l计算各结点初始时刻水头值计算各结点初始时刻水头值l利用差分方程计算各结点利用差分方程计算各结点t1t1时刻水头值时刻水头值l利用边界条件计算边界结点水头值利用边界条件计算边界结点水头值l重复重复2 2、3 3步，直到计算出拟计算的各个时刻的水头值步，直到计算出拟计算的各个时刻的水头值显式格式（续显式格式（续3）本讲稿第二十六页，共四十九页算例（续算例（续4）在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200200米，时间步长为米，时间步长为0.250.25天，分别计算各节点天，分别计算各节点各各时刻的水头值。时刻的水头值。本讲稿第二十七页，共四十九页Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 m02010101010100.252012.5101010100.502013.7510.6251010100.752014.531 11.250 10.156 10101.002015.078 11.797 10.391 10.039 101.252015.488 12.266 10.654 10.117 10算例（续算例（续5）本讲稿第二十八页，共四十九页Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 m02010101010100.252020101010100.502010201010100.7520300 2010101.0020-10 101.252010算例（续算例（续6）如果如果 t=1,t=1,则则本讲稿第二十九页，共四十九页(6)算例：算例：隐式格式隐式格式在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200米，时间步长为米，时间步长为0.25天，用隐式差分格式计算各节点个时刻的水头天，用隐式差分格式计算各节点个时刻的水头值。值。本讲稿第三十页，共四十九页在这个例子中，在这个例子中，解：隐式格式一般方程为解：隐式格式一般方程为于是有于是有根据初始条件得根据初始条件得根据边界条件得根据边界条件得本讲稿第三十一页，共四十九页由初始条件和边界条件由初始条件和边界条件由此解得由此解得t1t1时刻的水头值为时刻的水头值为在上述方程中取在上述方程中取 n=0,n=0,可以得到计算可以得到计算t1t1时刻水头值的方程时刻水头值的方程所以上述方程变成所以上述方程变成同理，可计算同理，可计算t2t2时刻时刻的水头值的水头值本讲稿第三十二页，共四十九页2.地下水数值模型建模步骤地下水数值模型建模步骤本讲稿第三十三页，共四十九页模拟步骤模拟步骤o建立概念模型建立概念模型 o建立数学模型建立数学模型 o数值方法及软件（编数值方法及软件（编程）程）o参数参数 o拟合模拟：模型校正拟合模拟：模型校正与检验与检验 o参数敏感性分析参数敏感性分析 o预测模拟预测模拟软件软件本讲稿第三十四页，共四十九页一、概念模型（模型概化）一、概念模型（模型概化）根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、岩石矿物、根据详细的地形地貌、地质、水文地质、构造地质、水文地球化学、岩石矿物、水文、气象、工农业利用情况等水文、气象、工农业利用情况等模拟的区域模拟的区域:含水层类型含水层类型:潜水（无压）、承压、混合、多层潜水（无压）、承压、混合、多层维数维数:一维、二维、三维一维、二维、三维水流状态：稳定流水流状态：稳定流/非稳定流、饱和流非稳定流、饱和流/非饱和流非饱和流介质状况介质状况:均质和非均质均质和非均质/各向同性和各向异性各向同性和各向异性孔隙孔隙/裂隙裂隙/双重介质双重介质流体的密度差流体的密度差边界条件和初始条件边界条件和初始条件 必要时需进行一系列的室内试验与野外试验必要时需进行一系列的室内试验与野外试验,以获取有关参数以获取有关参数,如渗透系数、弥如渗透系数、弥散系数、分配系数、反应速率常数等。散系数、分配系数、反应速率常数等。本讲稿第三十五页，共四十九页收集研究区已有收集研究区已有水文地质资料水文地质资料概念模型概念模型边界概化边界概化内部结构概化内部结构概化完成模型概化图完成模型概化图地表水体地表水体水头边界水头边界流量边界流量边界断层接触边界断层接触边界隔水边界隔水边界水头边界水头边界流量边界流量边界岩体岩层接触边界岩体岩层接触边界天然分水岭天然分水岭据含水层组类型、据含水层组类型、结构、岩性结构、岩性含水介质含水介质地下水运动状态地下水运动状态水文地质参数水文地质参数确定层组的均质、非均质，确定层组的均质、非均质，各向同行、异性，稳定、非各向同行、异性，稳定、非稳定，潜水、承压水稳定，潜水、承压水空间分布空间分布导水、储水系数、主渗透方向导水、储水系数、主渗透方向目的层与相邻层关系目的层与相邻层关系均质、非均质，各向同、异性均质、非均质，各向同、异性层、紊流，二、三维层、紊流，二、三维时间概化时间概化渗透、储水系数，给水度及单位涌渗透、储水系数，给水度及单位涌水量，含水层分布规律，地下水流水量，含水层分布规律，地下水流场，水化学场、温度场的空间概化场，水化学场、温度场的空间概化平面图平面图剖面图剖面图研究范围、主要居民点、标志性地形、地貌、地表河流、研究范围、主要居民点、标志性地形、地貌、地表河流、湖泊、开采井、地下水天然露头、含水层的各类钻孔湖泊、开采井、地下水天然露头、含水层的各类钻孔地表地理要素、含水层结构、地质结构、地下水水位、地表地理要素、含水层结构、地质结构、地下水水位、各类源汇项及其性质各类源汇项及其性质本讲稿第三十六页，共四十九页二、数学模型二、数学模型 三维地下水流动问题控制方程三维地下水流动问题控制方程第二类边界条件第二类边界条件第一类边界条件第一类边界条件初始条件初始条件本讲稿第三十七页，共四十九页 绝大部分数学模型是无法用解析法求解的，数值绝大部分数学模型是无法用解析法求解的，数值化就是将数学模型转化为可解的数值模型。化就是将数学模型转化为可解的数值模型。三、数值方法及软件（或编程）三、数值方法及软件（或编程）有限差分法有限差分法有限单元法有限单元法积分有限差分法积分有限差分法半解析半数值法半解析半数值法边界元法边界元法本讲稿第三十八页，共四十九页有限差分法有限差分法:MOFLOW系列系列GMS 中中MODFLOWVisual MODFLOWProcessing MODFLOW有限单元法有限单元法:FEFLOW积分有限差分法积分有限差分法:TOUGH2，TOUGH REACT软件软件本讲稿第三十九页，共四十九页四、模型参数四、模型参数o含水层参数：渗透系数，弹性释水系数（重力给水度）含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等，孔隙度等o源汇项：源汇项：n大气降水入渗系数（分区、数值）大气降水入渗系数（分区、数值）n蒸发排泄系数蒸发排泄系数n地表水体（河流、湖泊、水库等）水位、底面高程、底面岩地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等）性特征（厚度、渗透系数等）n渠系灌溉入渗系数渠系灌溉入渗系数n人工开采（点状、面状）人工开采（点状、面状）o边界条件边界条件o初始条件初始条件本讲稿第四十页，共四十九页参数的不确定性参数的不确定性钻孔太少，钻孔太少，地层资料少，地层资料少，钻孔多，含钻孔多，含水层结构会水层结构会发生变化发生变化本讲稿第四十一页，共四十九页五、模拟：模型校正（参数识别五、模拟：模型校正（参数识别）将模拟结果与实测结果比较将模拟结果与实测结果比较,进行参数调整进行参数调整,使模拟结使模拟结果在给定的误差范围内与实测结果吻合。果在给定的误差范围内与实测结果吻合。调参过程是一个复杂而辛苦的工作调参过程是一个复杂而辛苦的工作,所调整的参数必须所调整的参数必须符合模拟区的具体情况。符合模拟区的具体情况。人机交互与自动调参相结合。尽管自动调参程序人机交互与自动调参相结合。尽管自动调参程序(如如PEST PEST),),也不能代替人的工作。也不能代替人的工作。本讲稿第四十二页，共四十九页五、模拟：模型检验五、模拟：模型检验 模型验证是在模型校正的基础上模型验证是在模型校正的基础上,进一步进一步调整参数调整参数,使模拟结果与第二次实测结果吻合使模拟结果与第二次实测结果吻合,以进一步提高模型的置信度。以进一步提高模型的置信度。本讲稿第四十三页，共四十九页六、灵敏度分析六、灵敏度分析 校正后的模型受参数值的时空分布、边界校正后的模型受参数值的时空分布、边界条件、水流状态等条件、水流状态等不确定度不确定度的影响。灵敏度分的影响。灵敏度分析就是为了确定不确定度对校正模型的影响程析就是为了确定不确定度对校正模型的影响程度。度。本讲稿第四十四页，共四十九页七、预测七、预测 用校正的参数值进行预测用校正的参数值进行预测,预测时需估算未预测时需估算未来的水流状态。来的水流状态。本讲稿第四十五页，共四十九页后续检查与模型的重建（完善）后续检查与模型的重建（完善）后续检查在模拟研究结束数年后进行。收集新的野外数后续检查在模拟研究结束数年后进行。收集新的野外数据以确定预测结果是否正确。如果模拟结果精确据以确定预测结果是否正确。如果模拟结果精确,则该模则该模型对该模拟区来说是有效的。由于场址的唯一性型对该模拟区来说是有效的。由于场址的唯一性,故模故模型只对该模拟区有效。型只对该模拟区有效。后续检查应后续检查应在预测结束足够长的时间后进行在预测结束足够长的时间后进行,以便有以便有足够的时间发生明显的变化。足够的时间发生明显的变化。本讲稿第四十六页，共四十九页模型的再设计模型的再设计 一般来说一般来说,后续检查会发现系统性能的变后续检查会发现系统性能的变化化,从而导致概念模型和模型参数的修改。一从而导致概念模型和模型参数的修改。一般来说般来说,所有模拟研究都应该进行到第五步所有模拟研究都应该进行到第五步,即校正灵敏度分析。即校正灵敏度分析。本讲稿第四十七页，共四十九页地地裂裂缝缝自然地理及水文地质条件自然地理及水文地质条件边边界界敏敏感感性性分分析析数学模型数学模型网网格格化化概念模型概念模型介质类型、结构特征介质类型、结构特征地下水补、径、排特征地下水补、径、排特征数值模型数值模型数值模型的建立数值模型的建立边边界界条条件件初初始始条条件件源源汇汇项项降降雨雨入入渗渗地地表表水水体体入入渗渗灌灌溉溉、渠渠系系蒸蒸发发排排泄泄人人工工开开采采介介质质参参数数观观测测孔孔动动态态拟拟合合流流场场拟拟合合水水均均衡衡对对比比分分析析观观测测孔孔动动态态检检验验流流场场检检验验水水均均衡衡对对比比检检验验拟合调参拟合调参识识别别结结果果不不符符合合要要求求修修正正概概念念模模型型小结：地下水流动小结：地下水流动模型构建过程模型构建过程识识别别结结果果不不符符合合要要求求修修正正概概念念模模型型模型应用模型应用水水资资源源量量评评价价预预测测开开采采方方案案研研究究环环境境地地质质问问题题地地面面沉沉降降地地下下水水污污染染研研究究大大型型工工程程对对环环境境的的影影响响模型识别模型识别模型检验模型检验本讲稿第四十八页，共四十九页3.地下水流数值模拟资料需求地下水流数值模拟资料需求o空间展布相关资料空间展布相关资料n地表高程、分层数据地表高程、分层数据n边界位置边界位置o含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等含水层参数：渗透系数，弹性释水系数（重力给水度），孔隙度等o源汇项：源汇项：n大气降水入渗系数（分区、数值）大气降水入渗系数（分区、数值）n蒸发排泄系数蒸发排泄系数n地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等）地表水体（河流、湖泊、水库等）水位、底面高程、底面岩性特征（厚度、渗透系数等）n渠系灌溉入渗系数渠系灌溉入渗系数n人工开采（点状、面状）：开采井位置、井结构、开采量动态人工开采（点状、面状）：开采井位置、井结构、开采量动态o边界条件：边界类型、水头或流量边界条件：边界类型、水头或流量o初始条件：统测水位初始条件：统测水位o水位动态观测资料：观测孔位置、结构、水位时间变化水位动态观测资料：观测孔位置、结构、水位时间变化本讲稿第四十九页，共四十九页