第1-3章 基本原理

1、水文地质数值计算水文地质数值计算郭巧娜郭巧娜地球科学与工程学院地球科学与工程学院第一章第一章 地下水流动定解问题概述地下水流动定解问题概述1、三维流微分方程、三维流微分方程2、数学模型、数学模型 3、微分方程、微分方程4、边界条件、边界条件5、模型概化、模型概化6、实例、实例tzyvxtzyxx),(| )(方向流入tzyvxtzyxxx),(| )(方向流出),(| )(tzyxxxv),(| )(tzyxxv 取右图所示的微小六面取右图所示的微小六面体。设与体。设与x,y,z,方向对应的方向对应的主渗透系数分别为主渗透系数分别为Kxx, Kyy,Kzz；建立均衡期；建立均衡期 t时段内，微

2、小均衡六面体时段内，微小均衡六面体的水量守恒方程。的水量守恒方程。1、三维流微分方程、三维流微分方程1、三维流微分方程、三维流微分方程某一点为 , yyy , 某一点为xxx某一点为 , zzzx，y，z方向流入方向流入流出分别为流出分别为:l t时段内，六面体水量变化量为：时段内，六面体水量变化量为：tzyxzvyvxvzyx|)(|)(|)(tzyxxvtzyvvxxxxxx|tzxyyvtzxvvyyyyyy|tyxzzvtyxvvzzzzzz|六面体内地下水储存量的变化为六面体内地下水储存量的变化为由水均衡原理得由水均衡原理得1、三维流微分方程、三维流微分方程HzyxstHwzvyvx

3、vszyx)()()(Hzyxtzyxzvyvxvszyx|)(|)(|)(方程两端除以方程两端除以t，并取，并取,和和，则，则zHKzyHKyxHKxzzyyxxzHKvyHKvxHKvzzzyyyxxx一般密度的空间变化率很小，故一般密度的空间变化率很小，故于是有于是有()()()yxzvvvxyz由达西定律由达西定律有有（2）水流连续性方程左端项水流连续性方程左端项zHKzyHKyxHKxzzyyxx,都很小，可以忽略。都很小，可以忽略。1、三维流微分方程、三维流微分方程()esnM 其中为比储水系数)7(tHwzHKzyHKyxHKxszzyyxx上式为非均质各向异性承压含水层的偏微分

4、方程。上式为非均质各向异性承压含水层的偏微分方程。均质各向异性非稳定流均质各向异性非稳定流)10(222222tHwzHKyHKxHKszzyyxx0222222wzHKyHKxHKzzyyxx均质各向异性稳定流均质各向异性稳定流得到地下水三维流动微分方程得到地下水三维流动微分方程1/L1、三维流微分方程、三维流微分方程微分方程微分方程定解条件定解条件边界条件边界条件初始条件初始条件 已知已知t=0时的因变量，时的因变量，H(x,y,z,0)=H0(x,y,z)已知水头边界已知水头边界(I(I类边界类边界) )H(x,y,z,t)=f(x,y,z,t) (x,y,z) B1特例：特例：定水头边

5、界定水头边界 H(x,y,z,t)=C已知流量边界已知流量边界特例：隔水边界特例：隔水边界2B),(),(zyxtzyxfnH0nH2、数学模型、数学模型 数学模型数学模型地下水运动的数学模型结构地下水运动的数学模型结构(1) 第一类边界条件(简记为Dirichlet BC)n第一类边界条件也称为给定水头边界条件，指的是渗流区某一部分边界上，各点水头在某一时刻是已知的，表示为：n ),(),(1tzyxtzyxHDDtzyx),(2)第二类边界条件(简记为Neumann BC)n第二类边界条件也称为定流量边界条件，是指某一部分边界单位面积上流入(流出时值为负)的流量已知。相应的边界条件表示为：

6、),(1tzyxqnHKNNtzyx),(注水井或抽水井也可以作为内边界来处理，此时的承压含水层可忽略其水头的垂向变化，注水井或抽水井开口于整个厚度，如图所示，设井的半径为 ，单位时间抽水量为 ，井壁法向 和 轴方向相反，如图中所示。R)(tQ承压抽水井示意图承压抽水井示意图则由Darcy 定律，有：)()(2tQtKMRvrHKnHKt)(vRMtQrHKRr2)( （3） 第三类边界条件(Rolin BC)n第三类边界条件也称为混合边界条件，是指知道某段边界上 和 的线性组合：HnHn临海含水层与海水之间有一层薄淤泥层，则淤泥层两侧同一位置上的点存在水头差，内侧含水层水头 、渗透系数分别为

7、 ，淤泥层厚度为 ，水头、渗透系数分别为 和 ，外侧海水水头为 ，忽略淤泥层内孔隙水的弹性贮存的变化，在淤泥层与含水层交界面上根据流量的相等关系，有：HKmhKseaHn淤泥层与含水层交界面上根据流量的相等关系，mHHKnhKnHKseaseaHKmKHKmKnHKmKseaHKmK其中其中 （4）等值面边界条件n 等值面边界条件适用于不完全井的情况，即井壁和底部都有含水层中的水流渗入。等值面边界条件的表达形式是： 已知，并有正、负之分，分别对应于单位时间内的注入水量和抽出水量。 是井壁或渠壁与含水层的所有接触面。)(tQdsnHK )(tQ（5） 自由面边界条件n自由面边界条件即潜水面边界条

8、件，若潜水含水层中水头为 ，自由面方程： ，那么，在自由面上，有：),(tzyxH),(tyxz),(),(,(tyxttyxyxH（5） 自由面边界条件0),()tyxzezzzyyyyxxxxtnHKHKHK3、微分方程、微分方程tHwxHTxtHwyHTyxHTxtHwzHKzxHKxtHwzHKzyHKyxHKxxxyyxxszzxxszzyyxx一维流：平面二维流：剖面二维流：三维流： 数学模型数学模型4、边界条件、边界条件 边界条件边界条件：渗流区边界上水力特征，即边界上的水头分布和变：渗流区边界上水力特征，即边界上的水头分布和变化特征或流入流出含水层的水量分布和变化情况。主要有两

9、类：化特征或流入流出含水层的水量分布和变化情况。主要有两类： 1 已知边界上的水头分布规律已知边界上的水头分布规律 HB1=(x,y,t),其中其中(x,y,t)为已知函数。主要常见的是渗流为已知函数。主要常见的是渗流区与地表水体相接触。区与地表水体相接触。2已知边界上的单位宽度流量已知边界上的单位宽度流量q随时间的变化规律随时间的变化规律),(2tyxnHTB如已知流量为如已知流量为Q的承压含水层中完整的抽水井，其井壁可以看作的承压含水层中完整的抽水井，其井壁可以看作此类边界。此类边界。rQnHTB2 数学模型数学模型5例子例子:河间地块承压水流模型河间地块承压水流模型 设两条河流平行、完全

10、切割乘压含水层，含水层等厚、均质各向设两条河流平行、完全切割乘压含水层，含水层等厚、均质各向同性同性, ,无垂向补给或排泄，对于如图所示的坐标系，已知某时刻的含无垂向补给或排泄，对于如图所示的坐标系，已知某时刻的含水层各处的水头为水层各处的水头为2020米，自该时刻后，河水位分别为如图所示的函米，自该时刻后，河水位分别为如图所示的函数。试根据条件作合理简化建立其数学模型。数。试根据条件作合理简化建立其数学模型。（1 1）模型概化）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。由所述水文地质条件，可以概化为一维承压水流问题。（2 2）建立坐标系（如图）建立坐标系（如图） 取取x-x-

11、轴原点位于左端河，右侧为正向，设两河流间距为轴原点位于左端河，右侧为正向，设两河流间距为L. L. 纵轴为水头。纵轴为水头。（3 3）数学模型）数学模型22eHHKMtx)(00 xHHt)(),(210tHtHLXxLx 00;0tLx0t河间地块承压水流模型（续）河间地块承压水流模型（续） n1 有限差分法原理有限差分法原理n2 导数的有限差分近似表示导数的有限差分近似表示n3 承压一维流动有限差分法承压一维流动有限差分法n4 承压二维不稳定流有限差分法承压二维不稳定流有限差分法n5 源汇项的处理及井孔水头校正源汇项的处理及井孔水头校正n6不规则边界问题不规则边界问题n7 矩形变格距网格差

12、分矩形变格距网格差分第二章第二章 有限差分法有限差分法1 有限差分法的基本原理有限差分法的基本原理将连续的问题离散后求解：将连续的问题离散后求解： 方法一方法一以地下水流基本微分以地下水流基本微分方程及其定解条件为基础，方程及其定解条件为基础， 在在渗流区剖分基础上，用差商代渗流区剖分基础上，用差商代替微商，将地下水流微分方程替微商，将地下水流微分方程的求解转化为差分方程（代数的求解转化为差分方程（代数方程）求解。方程）求解。 方法二方法二在渗流区剖分的基础在渗流区剖分的基础上，直接由达西定律和水均衡上，直接由达西定律和水均衡原理，建立各个均衡区的水均原理，建立各个均衡区的水均衡方程，即差分方

13、程。衡方程，即差分方程。矩形网格矩形网格多边形网格多边形网格1 有限差分法的基本原理有限差分法的基本原理网格划分的基本类型网格划分的基本类型n（1）先划格线，）先划格线，格点位于网格中心格点位于网格中心均衡网格均衡网格节点网格节点网格o（2）先规定格点位）先规定格点位置，再垂直平分两相置，再垂直平分两相邻结点的连线作格线，邻结点的连线作格线，形成的网格即为水均形成的网格即为水均衡区衡区基本概念基本概念为时间步长。通常称为等或不相等的时段将连续的时间分割成相称为空间步长。方向上的格距和方向上的格距nnttx,yyx将时间离散点和空间离散点联合组成的网格称为时空网格。将时间离散点和空间离散点联合组

14、成的网格称为时空网格。有限差分的基本原理：某点处水头函数的导数用该点和几个相邻点处水头值及其间有限差分的基本原理：某点处水头函数的导数用该点和几个相邻点处水头值及其间距近似表示。距近似表示。MODFLOW网格系统网格系统x000()()()f xxf xfxx 导导数数的的有有限限差差商商近近似似导数的定义导数的定义 当当非常小的时候，有非常小的时候，有 上式右端项即为上式右端项即为f(x)(x)在在x0 0处的差商。处的差商。 这样定义的差商很容易理解，但不知道用差商代替微商所产生的误差。下面利用泰勒公式导出差商及其误差。xxfxxfxfx)()(lim)(0000方法一：方法一：差商代替微

15、商差商代替微商2 导数的有限差分近似表示导数的有限差分近似表示2000( )()()()()2!ff xxf xfxxx 2000( )()()()()2!ff xxf xfxxx 000()()()()f xxf xfxOxx 已知泰勒公式 由A得： AB 由B 得： 000()()()()f xf xxfxOxx 称称 为为f(x)在在x0处的处的一阶前向差商，一阶前向差商， 为为截断误差截断误差。xxfxxf)()(00)( xO 称称 为为f(x)在在x0处的处的一阶后向差商，一阶后向差商， 为为截断误差截断误差。xxxfxf)()(00)( xO 方法一方法一200002()2 ()

16、()()()()f xxf xf xxfxOxx 由A-B可以得： 由A+B可以得：2000()()()()2f xxf xxfxOxx (4)23400000()()( )()()()()()()2!3!4!fxfxff xxf xfxxxxx (4)23400000()()( )()()()()()()2!3!4!fxfxff xxf xfxxxxx AB称称 为为f(x)在在x0处的处的一阶中心差商，一阶中心差商， 为为截断误差截断误差。xxxfxxf2)()(002)( xO 称称 为为f(x)在在x0处的处的二阶二阶中心中心差商，差商， 为为截断误差截断误差。2000)()()(2)

17、(xxxfxfxxf2)( xO 方法一方法一n对于偏导数（偏微商），类似可以得到相应的差商对于偏导数（偏微商），类似可以得到相应的差商：ttxHttxHttxH),(),(),(000000 xtxHtxxHxtxH),(),(),(00000020000002002)(),(),(2),(),(xtxxHtxHtxxHxtxH方法一方法一),(|tzyxxxv),(|tzyxxvtzyvxtzyxx),(| )(方向流入tzyvxtzyxxx),(| )(方向流出 取右图所示得微小六面体。设与取右图所示得微小六面体。设与x,y,z,方向对应得主渗透系数分别为方向对应得主渗透系数分别为Kx,

18、 Ky,Kz；建立均衡期；建立均衡期 t时段内，时段内，微小均衡六面体的水量守恒方程。微小均衡六面体的水量守恒方程。方法二：达西定律和水均衡原理基于基于达西定律达西定律，x，y，z方向流入方向流入流出分别为流出分别为:l t时段内，侧向流入与源汇项导致六面体水量变化量为：时段内，侧向流入与源汇项导致六面体水量变化量为：tyxzyHHyHHtzxvvkjikjikjikjiyyyyy|, 1, 1,tyxzxHHxHHtzyvvkjikjikjikjixxxxx|, 1, 1tyxzzHHzHHtyxvvkjikjikjikjizzzzz|,1,1,ABC方法二：达西定律和水均衡原理DtyxzA

19、+B+C+D源汇项源汇项六面体内地下水储存量的变化为六面体内地下水储存量的变化为由水均衡原理得三维地下水流动方程的有限差分格式由水均衡原理得三维地下水流动方程的有限差分格式),(),(tzyxHttzyxHzyxHzyxssnkjinkjiskjikjikjikjikjikjikjikjikjikjikjikjiHHzyxtyxztyxzzHHzHHtyxzyHHyHHtyxzxHHxHH,1,1,1, 1, 1, 1, 1方法二：达西定律和水均衡原理有限差分法:三维（MODFLOW）差商代替微商差商代替微商3 承压一维流动有限差分法ttxhxtxhT),(),(22ttxHttxHxtxxH

20、txHtxxH),(),()(),(),(2),(00002000000ttxHttxHxttxxHttxHttxxH),(),()(),(),(2),(00002000000显式差分格式显式差分格式隐式差分格式隐式差分格式方法一方法一1, 3 , 2 )(21211 nxitHHxHHHTninininini控制方程控制方程1, 3 , 2 )(21211111 nxitHHxHHHTninininini网格剖分网格剖分nx个个算例：显式有限差格式 设两条河流平行、完全切割含水层，含水层等厚、均质设两条河流平行、完全切割含水层，含水层等厚、均质各向同性。各向同性。应用实例：河间地块承压水流模

21、型应用实例：河间地块承压水流模型3 承压一维流动有限差分法承压一维流动有限差分法步骤：（1 1）基础资料的分析）基础资料的分析（2 2）概念模型）概念模型（3 3）数学模型）数学模型（4 4）数值方法及计算机程序）数值方法及计算机程序（5 5）参数）参数（6 6）结果分析）结果分析 3 承压一维流动有限差分法承压一维流动有限差分法建立数学模型0,0tLxLx 022eHHKMtx)(00 xHHt)(),(210tHtHLXx（1 1）模型概化）模型概化 由所述水文地质条件，可以概化为一维承压水流问题。由所述水文地质条件，可以概化为一维承压水流问题。（2 2）建立坐标系）建立坐标系（如图），将

22、地下水流动系统空间结构放在坐标系内，（如图），将地下水流动系统空间结构放在坐标系内，从而量化各变量的取值范围。本例，取从而量化各变量的取值范围。本例，取x-x-轴原点位于左端河，右侧轴原点位于左端河，右侧为正向，设两河流间距为为正向，设两河流间距为L. L. （3 3）数学模型）数学模型0t3 承压一维流动有限差分法承压一维流动有限差分法差分方程及其解法显式格式将（0L）分成 N 等份，Nlx 1)1)网格剖分：网格剖分：xixittntn取时间步长取时间步长 ，记，记 （n n=0=0、1 1、2 2、3 3、44）记 ，(i=0,1,2,3，4N)2 2）建立差分方程）建立差分方程：在网格

23、系统中任意取一点：在网格系统中任意取一点),(nitx),(txH),(nitx设设是问题的解，则在是问题的解，则在处有处有记为（记为（i i，n n）22( , )nneiniiHHKMx ttx3 承压一维流动有限差分法承压一维流动有限差分法22( , )nneiniiHHKMx ttx1()nnniiiHHHOttt2211222()()nnnniiiiHHHHOxxx用差商代替微商用差商代替微商： 将上述两式舍去余项，代入方程并记将上述两式舍去余项，代入方程并记niHnih为为11122(,)()nnnnniiiiiinehhhhhKMx txt2)()(xOtO显然该式具有截断误差得

24、到得到显式格式（续1）111(12 )( , )nnnniiiiinethhhhx t引入无量纲变量引入无量纲变量：将该式子代入得到：将该式子代入得到： （i=1，2，3，.N-1），（n=1，2，3，.） 显式格式（续2）2)( xtKMe11122(,)()nnnnniiiiiinehhhhhKMx txt3）显示差分方程的求解n计算各结点初始时刻水头值计算各结点初始时刻水头值n利用差分方程计算各结点利用差分方程计算各结点t1t1时刻水头值时刻水头值n利用边界条件计算边界结点水头值利用边界条件计算边界结点水头值n重复重复2 2、3 3步，直到计算出拟计算的各个时刻的水头步，直到计算出拟计算

25、的各个时刻的水头值值)(00 xHHt显式格式（续显式格式（续3）22eHHKMtx)(00 xHHt)(),(210tHtHLXx算例（续4）004. 0em10)(m,20)( m,20 m/d,821ttMK在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200200米，时间步长为米，时间步长为0.250.25天，分别计算各节点天，分别计算各节点各各时刻的时刻的水头值。水头值。0 x20m,0 x10m,)(0 xH4/1)200(*004. 025. 0*20*8)(22xtKMeTime/dayx=0 mx=200 mx=400 mx=600 mx=800 m

26、x=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 10nieninininithhhh,111/)21 ()2(4141)4121 (4111111nininininininihhhhhhh算例（续5）Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 m020101010101

27、00.252020101010100.502010201010100.7520300 2010101.0020-10 101.252010算例（续6）1)200(*004. 01*20*8)(22xtKMenieninininithhhh,111/)21 (nininininininihhhhhhh11111) 121 (如果如果 t=1,t=1,则则算例：隐式格式22eHHKMtx)(00 xHHt)(),(210tHtHLXx004. 0e在上述模型中，设在上述模型中，设L=1000米米取空间步长为取空间步长为200米，时间步长为米，时间步长为0.25天，用隐式差分格式计算各节天，用隐式差

28、分格式计算各节点个时刻的水头值。点个时刻的水头值。m10)(m,20)( 10m,)( m,20 m/d,821ttxHMK11122332210111112131211/./.2121.212121nNnNenNnNenNnennennnennNnNnnnhththththhthhhhhh4/1)200(004. 025. 0208)(22xtKMe在这个例子中，在这个例子中，.0,1,2,3,4, m,10,m2050nhhnnm10,m10,m10,m1004030201hhhh154321011413121125. 025. 0225. 00025. 0225. 00025. 0225

29、. 00025. 02nnnnnnnnnnhhhhhhhhhh解：隐式格式一般方程为解：隐式格式一般方程为于是有于是有根据初始条件得根据初始条件得根据边界条件得根据边界条件得.0,1,2,3,4, m,10,m2050nhhnnm10,m10,m10,m1004030201hhhh1504030210011413121125. 025. 0225. 00025. 0225. 00025. 0225. 00025. 02hhhhhhhhhh由初始条件和边界条件由初始条件和边界条件5 .121010151025. 01010102025. 010225. 00025. 0225. 00025. 0

30、225. 00025. 0214131211hhhh由此解得由此解得t1t1时刻的水头值为时刻的水头值为m10,m10,m10,m1014131211hhhh在上述方程中取在上述方程中取 n=0,n=0,可以得到计算可以得到计算t1t1时刻水头值的方程时刻水头值的方程所以上述方程变成所以上述方程变成同理，可计算同理，可计算t2t2时刻的时刻的水头值水头值差分方程求解n一维显式差分格式 )(2)(2 )(2)(2111212212212331324322122321tHHxHHHtHHxHHHtHHxHHHtHHxHHHnninninninninninninninninninninnnnnnnn

31、nn网格个数为网格个数为ni直接求解直接求解差分方程求解n一维隐式差分格式一维隐式差分格式 )(2)(2 )(2)(2111211112212211121331321413122122131211tHHxHHHtHHxHHHtHHxHHHtHHxHHHnninninninninninninninninninninnnnnnnnnn网格个数为网格个数为ni迭代求解迭代求解方方程程组组PCGSIPSORWHSSAMGGMGMODFLOW差分方程的收敛性和稳定性n截断误差截断误差：用差商代替微商时，地下水流动方程产生用差商代替微商时，地下水流动方程产生的误差为截断误差。的误差为截断误差。n收敛性收敛

32、性：当空间步长和时间步长趋于当空间步长和时间步长趋于0 0时，有限差分方时，有限差分方程的精确解趋于地下水流动问题微分方程定解问题的程的精确解趋于地下水流动问题微分方程定解问题的精确解。则称该差分格式是收敛的。精确解。则称该差分格式是收敛的。n稳定性稳定性：如果在求解差分方程过程中，某时间步引入如果在求解差分方程过程中，某时间步引入某个误差，而在以后的各时段计算中，该误差不再扩某个误差，而在以后的各时段计算中，该误差不再扩大，则称该差分格式是稳定的。大，则称该差分格式是稳定的。102一维一维显示格式显示格式的收敛条件和稳定条件是：的收敛条件和稳定条件是：2()()OtOx第二类边界条件的处理n

33、前面介绍的几种差分格式，都是以第一类边界条件为例，说明其具体计算方法的。如果是第二类边界条件，怎么使用这些格式进行计算呢？为此，我们以下面的定解问题来讨论这个问题。设问题数学模型为出量。入量；取负值时，为流当流量取正值时，为流处的单宽流量。分别表示左、右两边界其中)(),(0)()(0)(0,0212100022tqtqttqxHTtqxHTLxxHHtLxtHxHTLxxt)()(),.,2 , 1 , 0()(,)(),.,2 , 1 , 0()()(0,211102121010100nnNnNnnnnnNnNnNnnNnnnnnnLxxtqTxhhtqTxhhMntqxhhTxHtqxH

34、TLxMntqxhhTxHtqxHTxxHxH得用后向差分公式代替，有对于右边界，得用前向差分公式代替处，有。对于左边界用后向差分公式代替用前向差分公式代替界条件中的偏导数，即用单侧差分公式代替边用一阶中心差分公式代替边界条件中的偏导数,.)1 , 0()(2)(2)(2)(2-,21111121111110111ntqTxhhtqTxhhtqxhhTtqxhhTxxxxxxxxnnNnNnnnnnNnNnnnNNN于是由边界条件得他们满足和为的边界点到虚线位置，记新扩充为此，首先将边界扩充4 承压二维不稳定流有限差分法n我们考虑无垂向补给、排泄、均质、各向同性、等厚的承压流动问题，边界条件第

35、一类，渗流区为矩形，则可用下述定解问题描述),(),(),(),()0 ;0(),()0 ;0(),(210210002222txHtxHtyHtyHbyaxyxHHbyaxtyxyHxHKMtHayyaxxte4 承压二维不稳定流有限差分法（1）显式格式-网格剖分.)3 , 2 , 1 , 0(,),.2 , 1 , 0()b, 0(),.2 , 1 , 0(), 0(ntnttNjyjyNbyNNixixNaxNanyjyyxixx记取时间步长等份，记分成将等份，记分成将（1）显式格式2）建立差分方程：22,1,21,1,2, 1, 1,22,22,)()()(),()(2)(2),(),

36、(),(),(),(yxtthhtyxyhhhxhhhKMtyxyHxHKMtHtxtxHnjityxnjinjienjinjinjinjinjinjinjinjinjinjinjieninji显然该式具有截断误差用差商代替微商，得到处有是问题的解，则在设记为点在网格系统中任意取一（1）显式格式.)3 , 2 , 1 , 0() 1N,.3 , 2 , 1, 1.N, 3 , 2 , 1(),()221 ()(,)(yx,1,1, 1, 11,22njityxthhhhhhytKMxtKMnjienjiyxnjiynjiynjixnjixnjieyex将该式子代入得到：引入无量纲变量：（1）显

37、式格式3）二维显示差分方程的求解计算各结点初始时刻水头值计算t1时刻水头值计算边界结点水头值重复2、3步，直到计算出拟计算的各个时刻水头值（1）显式格式n差分方程收敛性21)()(2yxxoto稳定条件是隐式格式的收敛条件和截断误差（2）隐式格式.)3 , 2 , 1; 1,.3 , 2 , 1; 1,.2 , 1(),()221 (-)(,)()()()(),()(2)(2),(1,1,1, 1, 122221 -,21,1,2, 1, 1,22,22,nNyjNxityxthhhhhhytKMxtKMyxtthhtyxyhhhxhhhKMtyxyHxHKMtHnjienjinjiyxnji

38、ynjiynjixnjixeyexnjinjienjinjinjinjinjinjinjinjinjinjinjie将该式子代入得到：引入无量纲变量：显然该式具有截断误差用差商代替微商，得到对于（2）隐式格式4）差分方程的收敛性和稳定性条件0, 0)()(2yxxoto稳定条件是隐式格式的收敛条件和截断误差5 源汇项的处理及井孔水头校正越流、入渗和抽水井等问题的处理越流、入渗和抽水井等问题的处理n如果考虑垂直渗流项（即源汇项），则二维承压流动如果考虑垂直渗流项（即源汇项），则二维承压流动微分方程可写成微分方程可写成n建立差分方程时，在结点处应加上这一项，它可具体建立差分方程时，在结点处应加上这

39、一项，它可具体表示为：表示为：2222HHHTTxy,()i juzi ji ji ji jzQKhhWMx y 为抽水或注水流量。；为入渗强度或蒸发强度为主含水层水头；为相邻含水层水头；为相邻弱透水层厚度；透系数；为相邻弱透水层垂向渗流等）；越流、入渗或蒸发、井为垂向渗流强度（包括jijijiujizzjiQWhhMK,井水位校正对比图井水位校正对比图 图图3-53-5（a a）初始网格）初始网格, ,()22wi j kwQxhhInTr有限差分法计算井水位的校正有限差分法计算井水位的校正 图图3-53-5（b b）加密网格）加密网格6不规则边界问题不规则边界问题l当研究区的几何形状属于简当研究区的几何形状属于简单形式（如矩形渗流区）时，单形式（如矩形渗流区）时，差分网格的划分往往将结点差分网格的划分往往将结点设在边界上。设在边界上。l然而，对于实际问题来说，然而，对于实际问题来说，边界通常不是那么规则，边边界通常不是那么规则，边界的某些部分，甚至大部分界的某些部分，甚至大部分不能与结点重合，我们称这不能与结点重合，我们称这种边界是不规则的。种边界是不规则的。l关于不规则边界问题，直接关于不规则边