地下水数值模拟03

1、第三章第三章 有限差分法有限差分法主要内容主要内容 有限差分法的基本原理有限差分法的基本原理 几种主要的差分格式几种主要的差分格式 二维渗流问题的差分方程二维渗流问题的差分方程 一般差分方程组的解法一般差分方程组的解法第一节第一节有限差分法的有限差分法的基本原理基本原理一、有限差分法的基本思想一、有限差分法的基本思想 1、基本原理、基本原理 从物理现象引出相应从物理现象引出相应微分方程微分方程（方程（方程+边界条件）边界条件）； 用差分网格用差分网格离散离散求解域；求解域； 用差分公式将基本方程转化为用差分公式将基本方程转化为差分方程差分方程（代数方程）（代数方程）； 用差分方程的解作为微分方

2、程的近似解。用差分方程的解作为微分方程的近似解。一、有限差分法的基本思想一、有限差分法的基本思想 2、求解步骤、求解步骤选取网格；选取网格；对微分方程及定解条件选择差分近似，列出对微分方程及定解条件选择差分近似，列出差分格式；差分格式；求解差分格式；求解差分格式；讨论差分格式解对于微分方程解的收敛性及讨论差分格式解对于微分方程解的收敛性及误差估计。误差估计。差分网格差分网格xyx，y空间步长空间步长t时间步长时间步长结点格点一、有限差分法的基本思想一、有限差分法的基本思想 3、优缺点、优缺点优优 点点缺缺 点点 1. 1.简单问题的数学表达简单问题的数学表达式和计算的执行过程式和计算的执行过程

3、比较直观、易懂；比较直观、易懂； 2.2.算法效率比较高，易算法效率比较高，易编程；编程；1 1、对自然边界处理的灵、对自然边界处理的灵活性较差。活性较差。2 2、对溶质运移等问题，、对溶质运移等问题，精度受限精度受限二、导数的有限差分近似表示二、导数的有限差分近似表示 1、差分的概念、差分的概念xxfxxfxydxdyxx)()(limlim00设有x的解析函数y(x)，函数y对x的导数为： 是函数对自变量的导数，又称微商。是函数对自变量的导数，又称微商。dxdy 、 分别称为函数及其自变量的差分分别称为函数及其自变量的差分yxdy、dx分别是函数及自变量的微分分别是函数及自变量的微分,xy

4、 为函数对自变量的差商。为函数对自变量的差商。由导数（微商）和差商的定义可知，由导数（微商）和差商的定义可知，当自变量的差分（增量）趋近于零时，当自变量的差分（增量）趋近于零时，就可以由差商得到导数。因此在数值就可以由差商得到导数。因此在数值计算中常用差商近似代替导数。计算中常用差商近似代替导数。二、导数的有限差分近似表示二、导数的有限差分近似表示 2、导数的有限差分形式、导数的有限差分形式用泰勒级数展开可以推导出导数的有限差分形式。用泰勒级数展开可以推导出导数的有限差分形式。333222! 3)(! 2)()()(dxfdxdxfdxdxdfxxfxxf333222! 3)(! 2)()()

5、(dxfdxdxfdxdxdfxxfxxf差分公式对比名称名称公式公式截断误差截断误差一一阶阶导导数数前差前差后差后差中心差中心差二阶导数二阶导数xxfxxfdxdf)()(xxxfxfdxdf)()(xxxfxxfdxdf2)()()(2xO )( xO )( xO 222)()()(2)(xxxfxfxxfdxfd)(2xO 思考与练习 高阶导数的有限差分形式高阶导数的有限差分形式33dxfd?yxf2?44xf? dxxdfdxxxdfdxxxdfxxxfxxfxxfdxddxfddxddxfd21 2222233 332222222222122221xxxfxxfxxfxxfxxfxx

6、fxxfxxfxxxxfxxfxxxfxfxxfxxfx思考与练习 高阶导数的有限差分形式高阶导数的有限差分形式yxf2?867527865310241 2221 2ffffhhffhffhhyfyfyxf015234678910111211931042031231102109220223221220444461 22221 2fffffhhfffhfffhfffhhxfxfxfxf三、简单水文地质模型的有限差分方程组三、简单水文地质模型的有限差分方程组 1、水文地质模型、水文地质模型u 以一维河间地块承压含水层中的水流问题为例。以一维河间地块承压含水层中的水流问题为例。u 考虑一个以通过考虑

7、一个以通过x0和和xL处的长且直的河流为界的承压含水层，该含水处的长且直的河流为界的承压含水层，该含水层均质各向同性，顶底板水平，上覆弱透水层，垂向补给强度为层均质各向同性，顶底板水平，上覆弱透水层，垂向补给强度为W，两河流，两河流边界的水位分别边界的水位分别0、 L为且不随时间变化。为且不随时间变化。u 试研究含水层的水头分布。试研究含水层的水头分布。 三、简单水文地质模型的有限差分方程组三、简单水文地质模型的有限差分方程组 2、数学模型、数学模型LLxxxHxHLxWxHT)()()0( 00022解析解212222CxCxTWHCxTWxHTWxHLLCLTWLx0122时，当020Cx

8、时，当02012CLTWLCL 00222xTWLLxTWxHL三、简单水文地质模型的有限差分方程组三、简单水文地质模型的有限差分方程组 3、有限差分方程、有限差分方程LLxxxHxHLxWxHT)()()0( 00022（1）网格剖分）网格剖分 沿河流的方向取单宽0,L作为研究区域，将L等分为n份，空间步长 对每个结点进行编号，结点编号由左向右依次为0，1，2，i，n。共有n+1个结点，其中2个边界结点，n-1个内结点。nLx 三、简单水文地质模型的有限差分方程组三、简单水文地质模型的有限差分方程组 3、有限差分方程、有限差分方程（2）建立有限差分方程）建立有限差分方程先任取一结点i进行分析

9、。 022WxHT0)(2211WxhhhTiii2112xTWhhhiii移项整理，得：移项整理，得： 对于所有内结点内结点1、2、n-1，建立结点相应的差分方程组 LnnnnnTxWhhTxWhhhTxWhhhTxWhh2122123232102212 2 2 2n-1个线性代数方程，未知量共n-1个故方程可解。 三、简单水文地质模型的有限差分方程组三、简单水文地质模型的有限差分方程组 3、有限差分方程、有限差分方程（3）求解）求解LnnnnnTxWhhTxWhhhTxWhhhTxWhh2122123232102212 2 2 2LnnTxWTxWTxWTxWhhhh/2112112112

10、222021221将方程表示成矩阵形式，则有： 线性代数方程组！线性代数方程组！ 系数矩阵中的元素仅位于三条对角线上，系数矩阵对称且系数矩阵中的元素仅位于三条对角线上，系数矩阵对称且正定，故称为正定，故称为三对角线性代数方程组。三对角线性代数方程组。 最有效的求最有效的求解方法解方法追赶法追赶法第二节第二节几种主要的差分格式几种主要的差分格式水文地质模型描述水文地质模型描述u 以一维河间地块承压含水层中的水流问题为例。u 含水层均质各向同性，不考虑垂向补给，两河流边界的水位随时间变化，分别0(t)、 L(t) 。u 试研究含水层的水头分布。 )0( )(),()0( )(),()0( )(),

11、()0 ,0( 0000*22TtttxHTtttxHLxxHtxHTtLxtHxHTLLxxt该问题属于一维承压非稳定流的定解问题。求解该问题，需要对空间该问题属于一维承压非稳定流的定解问题。求解该问题，需要对空间和时间进行离散，形成的差分网格称为时空网格。和时间进行离散，形成的差分网格称为时空网格。 网格剖分网格剖分以等距剖分为例以等距剖分为例将研究区域将研究区域0，L用直线等分为用直线等分为n份，把时间段份，把时间段0,T用直线等分成用直线等分成m份份以以 表示结点（表示结点（i,k）处的水头）处的水头 kiH)0( )(),()0( )(),()0( )(),()0 ,0( 0000*

12、22TtttxHTtttxHLxxHtxHTtLxtHxHTLLxxt 导数可以利用一阶、二阶导数的差商代替，由于一阶导数可以有三种差商表示，因此分别对水头关于时间的导数项分别运用前差、后差、中心差将得到三种差分格式。 显式有限差分显式有限差分 前差前差 隐式有限差分隐式有限差分 后差后差 中心式有限差分中心式有限差分 中心差中心差 一、一维显式有限差分格式一、一维显式有限差分格式tHxHT*22(i,k)(i-1,k)(i+1,k)(i,k+1)211)(2xhhhTkikikithhkiki1*向前差分向前差分=kikikikikihhhhhxtT1112*2)(整理得：整理得：2*)(

13、xtT定义定义kikikikihhhh111)21 (截断误差为：O(x2)+O(t) ktkhikktk+1hik+1，。 h11， h21 hn-11h01hl1,差分公（取k=1， i=1，i=2,i=n-1，t2t3 ，t4 ，t0h00， h10，hl0；差分方程，k=0，i=1，i=2,i=n-1,t1h11，hn-10。开始输入初始参数t=0对i=0,1,n循环执行xiHhi00t=t+t对i=1,n-1循环执行kikikikihhhh111)21 ( th010 thLl1输出结果：t，1ih10iihh tTsum?结束NoYes hik+1Hik+1?显式差分格式收显式差分

14、格式收敛和稳定的条件敛和稳定的条件2102)(Txt ，x, t：x, t：二、一维隐式有限差分格式二、一维隐式有限差分格式tHxHT*22(i,k+1)(i-1,k+1)(i+1,k+1)(i,k)211111)(2xhhhTkikikithhkiki1*向后差分向后差分=kikikikikihhhhhxtT1111112*2)(整理得：整理得：截断误差为：O(x2)+O(t) 2*)( xtT定义定义kikikikihhhh11111)21 (hik+1，。hik+1（i=1,i=2,i=l-1)差分方程代入形成由（l-1)个方程组成的线性方程组，联立求解。 112210111212111

15、2121212knknknkkknknkkhhhhhhhhhh开始输入初始参数t=0对i=0,1,n循环执行xiHhi00t=t+t对i=1,n-1循环线性方程组调用追赶法程序求解三对角线方程组输出结果：t，1ih10iihh tTsum?结束NoYeskikikikihhhh11111)21 ( 三、一维中心式有限差分格式三、一维中心式有限差分格式tk与与tk+1tk+t/2（i,k+1/2)(i,k+1/2)(i-1,k+1/2)(i+1,k+1/2)(i,k+1)k时刻k+1时刻2t(i,k)k+1/2时刻虚拟时刻变量虚拟时刻变量怎么取？怎么取？tHxHT*22截断误差为：O(x2)+O (t2 ) 21*kitHkikikixHxHTxHT2212221222121中心差分中心差分thhxhhhxhhhTkikikikikikikiki1*211211111)(2)(22整理得：整理得：kikikikikikihhhhhh1111111)22()22(2*)( xtT定义定义hik+1，thhxhhhxhhhkikikikikikikiki1211211111T)(21)(2）（ 小结显式差分隐式差分中心式差分差分方程公式向前差分向前差分向后差分向后差分中心