第1-3章基本原理

水文地质数值计算郭巧娜地球科学与工程学院第一章地下水流动定解问题概述1、三维流微分方程2、数学模型3、微分方程4、边界条件5、模型概化6、实例tzyvxtzyxx),,,(|)(方向流入tzyvxtzyxxx),,,(|)(方向流出),,,(|)(tzyxxxv),,,(|)(tzyxxv取右图所示的微小六面体。设与x,y,z,方向对应的主渗透系数分别为Kxx,Kyy,Kzz；建立均衡期t时段内，微小均衡六面体的水量守恒方程。1、三维流微分方程1、三维流微分方程某一点为],[yyy],[某一点为xxx某一点为],[zzzx，y，z方向流入—流出分别为:t时段内，六面体水量变化量为：tzyxzvyvxvzyx]|)(|)(|)([tzyxxvtzyvvxxxxxx|]||[tzxyyvtzxvvyyyyyy|]||[tyxzzvtyxvvzzzzzz|]||[六面体内地下水储存量的变化为由水均衡原理得1、三维流微分方程HzyxstHwzvyvxvszyx])()()([Hzyxtzyxzvyvxvszyx]|)(|)(|)([方程两端除以Δt，并取Δｘ→０,Δｙ→０,Δｚ→０和Δｔ→０，则zHKzyHKyxHKxzzyyxxzHKvyHKvxHKvzzzyyyxxx一般密度的空间变化率很小，故于是有()()()yxzvvvxyz由达西定律有（2）水流连续性方程左端项zHKzyHKyxHKxzzyyxx,,都很小，可以忽略。1、三维流微分方程()esnM其中为比储水系数)7(tHwzHKzyHKyxHKxszzyyxx上式为非均质各向异性承压含水层的偏微分方程。均质各向异性非稳定流)10(222222tHwzHKyHKxHKszzyyxx0222222wzHKyHKxHKzzyyxx均质各向异性稳定流得到地下水三维流动微分方程[1/L]1、三维流微分方程微分方程定解条件边界条件初始条件已知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已知流量边界特例：隔水边界2B),,(),,,(zyxtzyxfnH0nH2、数学模型数学模型地下水运动的数学模型结构(1)第一类边界条件(简记为DirichletBC)第一类边界条件也称为给定水头边界条件，指的是渗流区某一部分边界上，各点水头在某一时刻是已知的，表示为：),,,(),,,(1tzyxtzyxHDDtzyx),,,((2)第二类边界条件(简记为NeumannBC)第二类边界条件也称为定流量边界条件，是指某一部分边界单位面积上流入(流出时值为负)的流量已知。相应的边界条件表示为：),,,(1tzyxqnHKNNtzyx),,,(注水井或抽水井也可以作为内边界来处理，此时的承压含水层可忽略其水头的垂向变化，注水井或抽水井开口于整个厚度，如图所示，设井的半径为，单位时间抽水量为，井壁法向和轴方向相反，如图中所示。R)(tQ承压抽水井示意图则由Darcy定律，有：)()(2tQtKMRvrHKnHKt)(vRMtQrHKRr2)(（3）第三类边界条件(RolinBC)第三类边界条件也称为混合边界条件，是指知道某段边界上和的线性组合：HnH临海含水层与海水之间有一层薄淤泥层，则淤泥层两侧同一位置上的点存在水头差，内侧含水层水头、渗透系数分别为，淤泥层厚度为，水头、渗透系数分别为和，外侧海水水头为，忽略淤泥层内孔隙水的弹性贮存的变化，在淤泥层与含水层交界面上根据流量的相等关系，有：HKmhKseaH淤泥层与含水层交界面上根据流量的相等关系，mHHKnhKnHKseaseaHKmKHKmKnHKmKseaHKmK其中（4）等值面边界条件等值面边界条件适用于不完全井的情况，即井壁和底部都有含水层中的水流渗入。等值面边界条件的表达形式是：已知，并有正、负之分，分别对应于单位时间内的注入水量和抽出水量。是井壁或渠壁与含水层的所有接触面。)(tQdsnHK)(tQ（5）自由面边界条件自由面边界条件即潜水面边界条件，若潜水含水层中水头为，自由面方程：，那么，在自由面上，有：),,,(tzyxH),,(tyxz),,()),,,(,,(tyxttyxyxH（5）自由面边界条件0),,(])[tyxzezzzyyyyxxxxtnHKHKHK3、微分方程tHwxHTxtHwyHTyxHTxtHwzHKzxHKxtHwzHKzyHKyxHKxxxyyxxszzxxszzyyxx一维流：平面二维流：剖面二维流：三维流：数学模型4、边界条件边界条件：渗流区边界上水力特征，即边界上的水头分布和变化特征或流入流出含水层的水量分布和变化情况。主要有两类：1．已知边界上的水头分布规律HB1=ψ(x,y,t),其中ψ(x,y,t)为已知函数。主要常见的是渗流区与地表水体相接触。2．已知边界上的单位宽度流量q随时间的变化规律),,(2tyxnHTB如已知流量为Q的承压含水层中完整的抽水井，其井壁可以看作此类边界。rQnHTB2数学模型5例子:河间地块承压水流模型设两条河流平行、完全切割乘压含水层，含水层等厚、均质各向同性,无垂向补给或排泄，对于如图所示的坐标系，已知某时刻的含水层各处的水头为20米，自该时刻后，河水位分别为如图所示的函数。试根据条件作合理简化建立其数学模型。（1）模型概化由所述水文地质条件，可以概化为一维承压水流问题。（2）建立坐标系（如图）取x-轴原点位于左端河，右侧为正向，设两河流间距为L.纵轴为水头。（3）数学模型22eHHKMtx)(00xHHt)(),(210tHtHLXxLx00;0tLx0t河间地块承压水流模型（续）1有限差分法原理2导数的有限差分近似表示3承压一维流动有限差分法4承压二维不稳定流有限差分法5源汇项的处理及井孔水头校正6不规则边界问题7矩形变格距网格差分第二章有限差分法1有限差分法的基本原理将连续的问题离散后求解：方法一．以地下水流基本微分方程及其定解条件为基础，在渗流区剖分基础上，用差商代替微商，将地下水流微分方程的求解转化为差分方程（代数方程）求解。方法二．在渗流区剖分的基础上，直接由达西定律和水均衡原理，建立各个均衡区的水均衡方程，即差分方程。矩形网格多边形网格1有限差分法的基本原理——网格划分的基本类型（1）先划格线，格点位于网格中心均衡网格节点网格（2）先规定格点位置，再垂直平分两相邻结点的连线作格线，形成的网格即为水均衡区基本概念为时间步长。通常称为等或不相等的时段将连续的时间分割成相称为空间步长。方向上的格距和方向上的格距nnttx,yyx将时间离散点和空间离散点联合组成的网格称为时空网格。有限差分的基本原理：某点处水头函数的导数用该点和几个相邻点处水头值及其间距近似表示。MODFLOW网格系统x000()()()fxxfxfxx导数的有限差商近似导数的定义当非常小的时候，有上式右端项即为f(x)在x0处的差商。这样定义的差商很容易理解，但不知道用差商代替微商所产生的误差。下面利用泰勒公式导出差商及其误差。xxfxxfxfx)()(lim)(0000方法一：差商代替微商2导数的有限差分近似表示'''2000()()()()()2!ffxxfxfxxx'''2000()()()()()2!ffxxfxfxxx000()()()()fxxfxfxOxx已知泰勒公式①由A得：AB②由B得：000()()()()fxfxxfxOxx称为f(x)在x0处的一阶前向差商，为截断误差。xxfxxf)()(00)(xO称为f(x)在x0处的一阶后向差商，为截断误差。xxxfxf)()(00)(xO方法一''200002()2()()()()()fxxfxfxxfxOxx③由A-B可以得：④由A+B可以得：'2000()()()()2fxxfxxfxOxx'''''(4)'23400000()()()()()()()()()2!3!4!fxfxffxxfxfxxxxx'''''(4)'23400000()()()()()()()()()2!3!4!fxfxffxxfxfxxxxxAB称为f(x)在x0处的一阶中心差商，为截断误差。xxxfxxf2)()(002)(xO称为f(x)在x0处的二阶中心差商，为截断误差。2000)()()(2)(xxxfxfxxf2)(xO方法一对于偏导数（偏微商），类似可以得到相应的差商：ttxHttxHttxH),(),(),(000000xtxHtxxHxtxH),(),(),(00000020000002002)(),(),(2),(),(xtxxHtxHtxxHxtxH方法一),,,(|tzyxxxv),,,(|tzyxxvtzyvxtzyxx),,,(|)(方向流入tzyvxtzyxxx),,,(|)(方向流出取右图所示得微小六面体。设与x,y,z,方向对应得主渗透系数分别为Kx,Ky,Kz；建立均衡期t时段内，微小均衡六面体的水量守恒方程。方法二：达西定律和水均衡原理基于达西定律，x，y，z方向流入—流出分别为:t时段内，侧向流入与源汇项导致六面体水量变化量为：tyxzyHHyHHtzxvvkjikjikjikjiyyyyy][]||[,,,1,,,,1,tyxzxHHxHHtzyvvkjikjikjikjixxxxx][]||[,,,,1,,,,1tyxzzHHzHHtyxvvkjikjikjikjizzzzz][]||[,,1,,,,1,,ABC方法二：达西定律和水均衡原理DtyxzA+B+C+D源汇项六面体内地下水储存量的变化为由水均衡原理得三维地下水流动方程的有限差分格式),,,(),,,(tzyxHttzyxHzyxHzyxssnkjinkjiskjikjikjikjikjikjikjikjikjikjikjikjiHHzyxtyxztyxzzHHzHHtyxzyHHyHHtyxzxHHxHH,,1,,,,1,,,,1,,,,,1,,,,1,,,,,1,,,,1][][][方法二：达西定律和