第1-3章基本原理资料课件

水文地质数值计算水文地质数值计算郭巧娜郭巧娜地球科学与工程学院地球科学与工程学院第一章第一章 地下水流动定解问题概述地下水流动定解问题概述1、三维流微分方程、三维流微分方程2、数学模型、数学模型 3、微分方程、微分方程4、边界条件、边界条件5、模型概化、模型概化6、实例、实例 取右图所示的微小六面取右图所示的微小六面体。设与体。设与x,y,z,方向对应的方向对应的主渗透系数分别为主渗透系数分别为Kxx,Kyy,Kzz；建立均衡期；建立均衡期 t时段内，微小均衡六面体时段内，微小均衡六面体的水量守恒方程。的水量守恒方程。1、三维流微分方程、三维流微分方程1、三维流微分方程、三维流微分方程x，y，z方向流入方向流入流出分别为流出分别为:l t时段内，六面体水量变化量为：时段内，六面体水量变化量为：六面体内地下水储存量的变化为六面体内地下水储存量的变化为由水均衡原理得由水均衡原理得1、三维流微分方程、三维流微分方程方程两端除以方程两端除以t，并取，并取,和和，则，则一般密度的空间变化率很小，故一般密度的空间变化率很小，故于是有于是有由达西定律由达西定律有有（2）水流连续性方程左端项水流连续性方程左端项都很小，可以忽略。都很小，可以忽略。1、三维流微分方程、三维流微分方程上式为非均质各向异性承压含水层的偏微分方程。上式为非均质各向异性承压含水层的偏微分方程。均质各向异性非稳定流均质各向异性非稳定流均质各向异性稳定流均质各向异性稳定流得到地下水三维流动微分方程得到地下水三维流动微分方程1/L1、三维流微分方程、三维流微分方程微分方程微分方程定解条件定解条件边界条件边界条件初始条件初始条件已知已知t=0时的因变量，时的因变量，H(x,y,z,0)=H0(x,y,z)已知水头边界已知水头边界(I类边界类边界)H(x,y,z,t)=f(x,y,z,t)(x,y,z)B1特例：特例：定水头边界定水头边界 H(x,y,z,t)=C已知流量边界已知流量边界特例：隔水边界特例：隔水边界2、数学模型、数学模型 数学模型数学模型地下水运动的数学模型结构地下水运动的数学模型结构(1)第一类边界条件(简记为Dirichlet BC)n第一类边界条件也称为给定水头边界条件，指的是渗流区某一部分边界上，各点水头在某一时刻是已知的，表示为：n (2)第二类边界条件(简记为Neumann BC)n第二类边界条件也称为定流量边界条件，是指某一部分边界单位面积上流入(流出时值为负)的流量已知。相应的边界条件表示为：注水井或抽水井也可以作为内边界来处理，此时的承压含水层可忽略其水头的垂向变化，注水井或抽水井开口于整个厚度，如图所示，设井的半径为 ，单位时间抽水量为 ，井壁法向 和 轴方向相反，如图中所示。承压抽水井示意图承压抽水井示意图则由Darcy 定律，有：（3）第三类边界条件(Rolin BC)n第三类边界条件也称为混合边界条件，是指知道某段边界上 和 的线性组合：n临海含水层与海水之间有一层薄淤泥层，则淤泥层两侧同一位置上的点存在水头差，内侧含水层水头 、渗透系数分别为 ，淤泥层厚度为 ，水头、渗透系数分别为 和 ，外侧海水水头为 ，忽略淤泥层内孔隙水的弹性贮存的变化，在淤泥层与含水层交界面上根据流量的相等关系，有：n淤泥层与含水层交界面上根据流量的相等关系，其中其中（4）等值面边界条件n 等值面边界条件适用于不完全井的情况，即井壁和底部都有含水层中的水流渗入。等值面边界条件的表达形式是：已知，并有正、负之分，分别对应于单位时间内的注入水量和抽出水量。是井壁或渠壁与含水层的所有接触面。（5）自由面边界条件n自由面边界条件即潜水面边界条件，若潜水含水层中水头为 ，自由面方程：，那么，在自由面上，有：（5）自由面边界条件3、微分方程、微分方程 数学模型数学模型4、边界条件、边界条件 边界条件边界条件：渗流区边界上水力特征，即边界上的水头分布和变：渗流区边界上水力特征，即边界上的水头分布和变化特征或流入流出含水层的水量分布和变化情况。主要有两类：化特征或流入流出含水层的水量分布和变化情况。主要有两类：1 已知边界上的水头分布规律已知边界上的水头分布规律 HB1=(x,y,t),其中其中(x,y,t)为已知函数。主要常见的是渗流区为已知函数。主要常见的是渗流区与地表水体相接触。与地表水体相接触。2已知边界上的单位宽度流量已知边界上的单位宽度流量q随时间的变化规律随时间的变化规律如已知流量为如已知流量为Q的承压含水层中完整的抽水井，其井壁可以看作的承压含水层中完整的抽水井，其井壁可以看作此类边界。此类边界。数学模型数学模型5例子例子:河间地块承压水流模型河间地块承压水流模型 设两条河流平行、完全切割乘压含水层，含水层等厚、均质各向设两条河流平行、完全切割乘压含水层，含水层等厚、均质各向同性同性,无垂向补给或排泄，对于如图所示的坐标系，已知某时刻的含无垂向补给或排泄，对于如图所示的坐标系，已知某时刻的含水层各处的水头为水层各处的水头为2020米，自该时刻后，河水位分别为如图所示的函米，自该时刻后，河水位分别为如图所示的函数。试根据条件作合理简化建立其数学模型。数。试根据条件作合理简化建立其数学模型。（1 1）模型概化）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。由所述水文地质条件，可以概化为一维承压水流问题。（2 2）建立坐标系（如图）建立坐标系（如图）取取x-x-轴原点位于左端河，右侧为正向，设两河流间距为轴原点位于左端河，右侧为正向，设两河流间距为L.L.纵轴为水头。纵轴为水头。（3 3）数学模型）数学模型河间地块承压水流模型（续）河间地块承压水流模型（续）n1 有限差分法原理有限差分法原理n2 导数的有限差分近似表示导数的有限差分近似表示n3 承压一维流动有限差分法承压一维流动有限差分法n4 承压二维不稳定流有限差分法承压二维不稳定流有限差分法n5 源汇项的处理及井孔水头校正源汇项的处理及井孔水头校正n6不规则边界问题不规则边界问题n7 矩形变格距网格差分矩形变格距网格差分第二章第二章 有限差分法有限差分法1 有限差分法的基本原理有限差分法的基本原理将连续的问题离散后求解：将连续的问题离散后求解：方法一方法一以地下水流基本微分以地下水流基本微分方程及其定解条件为基础，方程及其定解条件为基础，在在渗流区剖分基础上，用差商代渗流区剖分基础上，用差商代替微商，将地下水流微分方程替微商，将地下水流微分方程的求解转化为差分方程（代数的求解转化为差分方程（代数方程）求解。方程）求解。方法二方法二在渗流区剖分的基础在渗流区剖分的基础上，直接由达西定律和水均衡上，直接由达西定律和水均衡原理，建立各个均衡区的水均原理，建立各个均衡区的水均衡方程，即差分方程。衡方程，即差分方程。矩形网格矩形网格多边形网格多边形网格1 有限差分法的基本原理有限差分法的基本原理网格划分的基本类型网格划分的基本类型n（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处的处的二阶二阶中心中心差商，差商，为为截断误差截断误差。方法一方法一n对于偏导数（偏微商），类似可以得到相应的差商：对于偏导数（偏微商），类似可以得到相应的差商：方法一方法一 取右图所示得微小六面体。设与取右图所示得微小六面体。设与x,y,z,方向对应得主渗透系数分别为方向对应得主渗透系数分别为Kx,Ky,Kz；建立均衡期；建立均衡期 t时段内，时段内，微小均衡六面体的水量守恒方程。微小均衡六面体的水量守恒方程。方法二：达西定律和水均衡原理基于基于达西定律达西定律，x，y，z方向流入方向流入流出分别为流出分别为:l t时段内，侧向流入与源汇项导致六面体水量变化量为：时段内，侧向流入与源汇项导致六面体水量变化量为：ABC方法二：达西定律和水均衡原理DA+B+C+D源汇项源汇项六面体内地下水储存量的变化为六面体内地下水储存量的变化为由水均衡原理得三维地下水流动方程的有限差分格式由水均衡原理得三维地下水流动方程的有限差分格式方法二：达西定律和水均衡原理有限差分法:三维（MODFLOW）差商代替微商差商代替微商3 承压一维流动有限差分法显式差分格式显式差分格式隐式差分格式隐式差分格式方法一方法一控制方程控制方程网格剖分网格剖分nx个个算例：显式有限差格式 设两条河流平行、完全切割含水层，含水层等厚、均质设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。各向同性。应用实例：河间地块承压水流模型应用实例：河间地块承压水流模型3 承压一维流动有限差分法承压一维流动有限差分法步骤：（1 1）基础资料的分析）基础资料的分析（2 2）概念模型）概念模型（3 3）数学模型）数学模型（4 4）数值方法及计算机程序）数值方法及计算机程序（5 5）参数）参数（6 6）结果分析）结果分析 3 承压一维流动有限差分法承压一维流动有限差分法建立数学模型（1 1）模型概化）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。由所述水文地质条件，可以概化为一维承压水流问题。（2 2）建立坐标系）建立坐标系（如图），将地下水流动系统空间结构放在坐标系内，（如图），将地下水流动系统空间结构放在坐标系内，从而量化各变量的取值范围。本例，取从而量化各变量的取值范围。本例，取x-x-轴原点位于左端河，右侧轴原点位于左端河，右侧为正向，设两河流间距为为正向，设两河流间距为L.L.（3 3）数学模型）数学模型3 承压一维流动有限差分法承压一维流动有限差分法差分方程及其解法显式格式将（0L）分成 N 等份，1)1)网格剖分：网格剖分：取时间步长取时间步长 ，记，记 （n n=0=0、1 1、2 2、3 3、4 4）记 ，(i=0,1,2,3，4N)2 2）建立差分方程）建立差分方程：在网格系统中任意取一点：在网格系统中任意取一点设设是问题的解，则在是问题的解，则在处有处有记为（记为（i i，n n）3 承压一维流动有限差分法承压一维流动有限差分法用差商代替微商用差商代替微商：将上述两式舍去余项，代入方程并记将上述两式舍去余项，代入方程并记为为显然该式具有截断误差得到得到显式格式（续1）引入无量纲变量引入无量纲变量：将该式子代入得到：将该式子代入得到：（i=1，2，3，.N-1），（n=1，2，3，.）显式格式（续2）3）显示差分方程的求解n计算各结点初始时刻水头值计算各结点初始时刻水头值n利用差分方程计算各结点利用差分方程计算各结点t1t1时刻水头值时刻水头值n利用边界条件计算边界结点水头值利用边界条件计算边界结点水头值n重复重复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,则则算例：隐式格式在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200米，时间步长为米，时间步长为0.25天，用隐式差分格式计算各节天，用隐式差分格式计算各节点个时刻的水头值。点个时刻的水头值。在这个例子中，在这个例子中，解：隐式格式一般方程为解：隐式格式一般方程为于是有于是有根据初始条件得根据初始条件得根据边界条件得根据边界条件得由初始条件和边界条件由初始条件和边界条件由此解得由此解得t1t1时刻的水头值为时刻的水头值为在上述方程中取在上述方程中取 n=0,n=0,可以得到计算可以得到计算t1t1时刻水头值的方程时刻水头值的方程所以上述方程变成所以上述方程变成同理，可计算同理，可计算t2t2时刻的时刻的水头值水头值差分方程求解n一维显式差分格式网格个数为网格个数为ni直接求解直接求解差分方程求解n一维隐式差分格式一维隐式差分格式网格个数为网格个数为ni迭代求解迭代求解方方程程组组PCGSIPSORWHSSAMGGMGMODFLOW差分方程的收敛性和稳定性n截断误差截断误差：用差商代替微商时，地下水流动方程产生用差商代替微商时，地下水流动方程产生的误差为截断误差。的误差为截断误差。n收敛性收敛性：当空间步长和时间步长趋于当空间步长和时间步长趋于0 0时，有限差分方时，有限差分方程的精确解趋于地下水流动问题微分方程定解问题的程的精确解趋于地下水流动问题微分方程定解问题的精确解。则称该差分格式是收敛的。精确解。则称该差分格式是收敛的。n稳定性稳定性：如果在求解差分方程过程中，某时间步引入如果在求解差分方程过程中，某时间步引入某个误差，而在以后的各时段计算中，该误差不再扩某个误差，而在以后的各时段计算中，该误差不再扩大，则称该差分格式是稳定的。大，则称该差分格式是稳定的。一维一维显示格式显示格式的收敛条件和稳定条件是：的收敛条件和稳定条件是：第二类边界条件的处理n前面介绍的几种差分格式，都是以第一类边界条件为例，说明其具体计算方法的。如果是第二类边界条件，怎么使用这些格式进行计算呢？为此，我们以下面的定解问题来讨论这个问题。设问题数学模型为用一阶中心差分公式代替边界条件中的偏导数4 承压二维不稳定流有限差分法n我们考虑无垂向补给、排泄、均质、各向同性、等厚的承压流动问题，边界条件第一类，渗流区为矩形，则可用下述定解问题描述4 承压二维不稳定流有限差分法（1）显式格式-网格剖分（1）显式格式2）建立差分方程：（1）显式格式（1）显式格式3）二维显示差分方程的求解计算各结点初始时刻水头值计算t1时刻水头值计算边界结点水头值重复2、3步，直到计算出拟计算的各个时刻水头值（1）显式格式n差分方程收敛性（2）隐式格式（2）隐式格式4）差分方程的收敛性和稳定性条件5 源汇项的处理及井孔水头校正越流、入渗和抽水井等问题的处理越流、入渗和抽水井等问题的处理n如果考虑垂直渗流项（即源汇项），则二维承压流动微如果考虑垂直渗流项（即源汇项），则二维承压流动微分方程可写成分方程可写成n建立差分方程时，在结点处应加上这一项，它可具体表建立差分方程时，在结点处应加上这一项，它可具体表示为：示为：井水位校正对比图井水位校正对比图 图3-53-5（a a）初始网格）初始网格有限差分法计算井水位的校正有限差分法计算井水位的校正 图3-53-5（b b）加密网格）加密网格6不规则边界问题不规则边界问题l当研究区的几何形状属于简当研究区的几何形状属于简单形式（如矩形渗流区）时，单形式（如矩形渗流区）时，差分网格的划分往往将结点差分网格的划分往往将结点设在边界上。设在边界上。l然而，对于实际问题来说，然而，对于实际问题来说，边界通常不是那么规则，边边界通常不是那么规则，边界的某些部分，甚至大部分界的某些部分，甚至大部分不能与结点重合，我们称这不能与结点重合，我们称这种边界是不规则的。种边界是不规则的。l关于不规则边界问题，直接关于不规则边界问题，直接取最靠近边界的曲折格线为取最靠近边界的曲折格线为近似边界近似边界河流、水库、湖泊等问题河流、水库、湖泊等问题7 矩形变格距网格差分n本节直接由达西定律和水均衡原理建立差分方程n达西定律：n水均衡原理：对某一研究对象，流入-流出=体系内质量变化量 研究对象可以是大区域，也可以是微分单元体。大区域的水均衡计算经常用于区域的水资源评价 表示当水头下降一个单位时，从单位体积饱水空隙介质中释放的水量或体积参考教材参考教材n周志芳主编，地下水运动数值模拟，自编教材，2003n薛禹群，谢春红，水文地质学的数值法，煤炭工业出版社，1980n陈崇希，唐仲华，地下水流动问题数值法，中陈崇希，唐仲华，地下水流动问题数值法，中国地质大学出版社，国地质大学出版社，1990n李俊亭，地下水数值模拟，地质出版社，1988n孙讷正，地下水流的数学模型和数值方法，地质出版社，1981n1 有限单元法原理有限单元法原理n2 承压二维不稳定流动问题的迦辽金方程承压二维不稳定流动问题的迦辽金方程n3 三角形单元迦辽金有限元法三角形单元迦辽金有限元法n4 矩形单元迦辽金有限元法矩形单元迦辽金有限元法n5 任意四边形等参有限元法任意四边形等参有限元法n6 无压流问题的有限元方法无压流问题的有限元方法n7 多层含水层越流系统准三维流问题多层含水层越流系统准三维流问题第三章第三章 迦辽金有限单元法迦辽金有限单元法1 有限单元法的基本原理有限单元法的基本原理有限元法的基本思想是将结构离散化，用有限个单元表示复杂对象，单元之间有限元法的基本思想是将结构离散化，用有限个单元表示复杂对象，单元之间通过节点相互连接，然后根据变形协调条件综合求解。通过节点相互连接，然后根据变形协调条件综合求解。有限单元法是求解偏微分方程定解问题的有效数值方法。有限单元法是求解偏微分方程定解问题的有效数值方法。与差分法类似，与差分法类似，有限元方法也是通过区域剖分和插值方法将描述地下水流动的有限元方法也是通过区域剖分和插值方法将描述地下水流动的定解问题化为代数方程组。定解问题化为代数方程组。依据建立代数方程组的途径不同，依据建立代数方程组的途径不同，有限单元法又分为迦辽金法和里兹法有限单元法又分为迦辽金法和里兹法。1 有限单元法的基本原理有限单元法的基本原理迦辽金有限单元法基本思想：迦辽金有限单元法基本思想：从剩余加权法出发，利用剖分插值的方法，从剩余加权法出发，利用剖分插值的方法，将数学模型的求解问题离散成常微分方程的将数学模型的求解问题离散成常微分方程的求解问题，在利用差分方法把常微分方程离求解问题，在利用差分方法把常微分方程离散成线性代数方程组的求解问题。散成线性代数方程组的求解问题。2 承压二维不稳定流动问题的迦辽金方程承压二维不稳定流动问题的迦辽金方程人有了知识，就会具备各种分析能力，人有了知识，就会具备各种分析能力，明辨是非的能力。明辨是非的能力。所以我们要勤恳读书，广泛阅读，所以我们要勤恳读书，广泛阅读，古人说古人说“书中自有黄金屋。书中自有黄金屋。”通过阅读科技书籍，我们能丰富知识，通过阅读科技书籍，我们能丰富知识，培养逻辑思维能力；培养逻辑思维能力；通过阅读文学作品，我们能提高文学鉴赏水平，通过阅读文学作品，我们能提高文学鉴赏水平，培养文学情趣；培养文学情趣；通过阅读报刊，我们能增长见识，扩大自己的知识面。通过阅读报刊，我们能增长见识，扩大自己的知识面。有许多书籍还能培养我们的道德情操，有许多书籍还能培养我们的道德情操，给我们巨大的精神力量，给我们巨大的精神力量，鼓舞我们前进鼓舞我们前进。