地下水数值模拟03-有限差分法

第三章有限差分法主要内容•有限差分法的基本原理•几种主要的差分格式•二维渗流问题的差分方程•一般差分方程组的解法第一节有限差分法的基本原理一、有限差分法的基本思想•1、基本原理•从物理现象引出相应微分方程（方程+边界条件）；•用差分网格离散求解域；•用差分公式将基本方程转化为差分方程（代数方程）；•用差分方程的解作为微分方程的近似解。一、有限差分法的基本思想•2、求解步骤选取网格；对微分方程及定解条件选择差分近似，列出差分格式；求解差分格式；讨论差分格式解对于微分方程解的收敛性及误差估计。差分网格xyΔx，Δy——空间步长Δt——时间步长结点格点一、有限差分法的基本思想•3、优缺点优点缺点1.简单问题的数学表达式和计算的执行过程比较直观、易懂；2.算法效率比较高，易编程；1、对自然边界处理的灵活性较差。2、对溶质运移等问题，精度受限二、导数的有限差分近似表示•1、差分的概念xxfxxfxydxdyxx)()(limlim00设有x的解析函数y(x)，函数y对x的导数为：——是函数对自变量的导数，又称微商。dxdy——、分别称为函数及其自变量的差分yx——dy、dx分别是函数及自变量的微分,xy——为函数对自变量的差商。由导数（微商）和差商的定义可知，当自变量的差分（增量）趋近于零时，就可以由差商得到导数。因此在数值计算中常用差商近似代替导数。二、导数的有限差分近似表示•2、导数的有限差分形式用泰勒级数展开可以推导出导数的有限差分形式。333222!3)(!2)()()(dxfdxdxfdxdxdfxxfxxf333222!3)(!2)()()(dxfdxdxfdxdxdfxxfxxf差分公式对比名称公式截断误差一阶导数前差后差中心差二阶导数xxfxxfdxdf)()(xxxfxfdxdf)()(xxxfxxfdxdf2)()()(2xO)(xO)(xO222)()()(2)(xxxfxfxxfdxfd)(2xO思考与练习•高阶导数的有限差分形式33dxfd?yxf2?44xf?dxxdfdxxxdfdxxxdfxxxfxxfxxfdxddxfddxddxfd212222233332222222222122221xxxfxxfxxfxxfxxfxxfxxfxxfxxxxfxxfxxxfxfxxfxxfx思考与练习•高阶导数的有限差分形式yxf2?86752786531024122212ffffhhffhffhhyfyfyxf015234678910111211931042031231102109220223221220444461222212fffffhhfffhfffhfffhhxfxfxfxf三、简单水文地质模型的有限差分方程组•1、水文地质模型以一维河间地块承压含水层中的水流问题为例。考虑一个以通过x＝0和x＝L处的长且直的河流为界的承压含水层，该含水层均质各向同性，顶底板水平，上覆弱透水层，垂向补给强度为W，两河流边界的水位分别φ0、φL为且不随时间变化。试研究含水层的水头分布。三、简单水文地质模型的有限差分方程组•2、数学模型LLxxxHxHLxWxHT)()()0(00022解析解212222CxCxTWHCxTWxHTWxHLLCLTWLx0122时，当020Cx时，当02012CLTWLCL00222xTWLLxTWxHL三、简单水文地质模型的有限差分方程组•3、有限差分方程LLxxxHxHLxWxHT)()()0(00022（1）网格剖分沿河流的方向取单宽[0,L]作为研究区域，将L等分为n份，空间步长对每个结点进行编号，结点编号由左向右依次为0，1，2，…，i，…，n。共有n+1个结点，其中2个边界结点，n-1个内结点。nLx三、简单水文地质模型的有限差分方程组•3、有限差分方程（2）建立有限差分方程先任取一结点i进行分析。022WxHT0)(2211WxhhhTiii2112xTWhhhiii移项整理，得：对于所有内结点1、2、…、n-1，建立结点相应的差分方程组LnnnnnTxWhhTxWhhhTxWhhhTxWhh2122123232102212222n-1个线性代数方程，未知量共n-1个故方程可解。三、简单水文地质模型的有限差分方程组•3、有限差分方程（3）求解LnnnnnTxWhhTxWhhhTxWhhhTxWhh2122123232102212222LnnTxWTxWTxWTxWhhhh////2112112112222021221将方程表示成矩阵形式，则有：线性代数方程组！系数矩阵中的元素仅位于三条对角线上，系数矩阵对称且正定，故称为三对角线性代数方程组。最有效的求解方法——追赶法第二节几种主要的差分格式水文地质模型描述以一维河间地块承压含水层中的水流问题为例。含水层均质各向同性，不考虑垂向补给，两河流边界的水位随时间变化，分别φ0(t)、φL(t)。试研究含水层的水头分布。)0()(),()0()(),()0()(),()0,0(0000*22TtttxHTtttxHLxxHtxHTtLxtHxHTLLxxt该问题属于一维承压非稳定流的定解问题。求解该问题，需要对空间和时间进行离散，形成的差分网格称为时空网格。网格剖分•以等距剖分为例•将研究区域[0，L]用直线等分为n份，把时间段[0,T]用直线等分成m份•以表示结点（i,k）处的水头kiH)0()(),()0()(),()0()(),()0,0(0000*22TtttxHTtttxHLxxHtxHTtLxtHxHTLLxxt导数可以利用一阶、二阶导数的差商代替，由于一阶导数可以有三种差商表示，因此分别对水头关于时间的导数项分别运用前差、后差、中心差将得到三种差分格式。显式有限差分←前差隐式有限差分←后差中心式有限差分←中心差一、一维显式有限差分格式tHxHT*22(i,k)(i-1,k)(i+1,k)(i,k+1)211)(2xhhhTkikikithhkiki1*向前差分=kikikikikihhhhhxtT1112*2)(整理得：2*)(xtT定义kikikikihhhh111)21(截断误差为：O([Δx]2)+O(Δt)只要知道了k时段初始时刻tk各结点的hik值，便可计算出k时段末了时刻tk+1的hik+1值，各方程可独立求解，因此，这种方程称为显式有限差分方程。差分方程与离散后的定界条件构成了数学模型的显式差分方程问题。其求解步骤如下：由以上计算的h11，h21…hn-11值及由边界条件计算的h01和hl1,再次利用差分公式（取k=1，i=1，i=2,…i=n-1，便可计算得t2时刻各结点的水头值。如此重复，便可计算出t3，t4，…各时刻的水头分布值。由初始条件给出t0时刻各结点的水头值h00，h10，…hl0；再根据差分方程，在k=0时，分别取i=1，i=2,…i=n-1,便可求得t1时刻各内结点的水头值h11，…hn-10。显式差分方程的求解开始输入初始参数t=0对i=0,1,…,n循环执行xiHhi00t=t+Δt对i=1,…,n-1循环执行kikikikihhhh111)21(th010thLl1输出结果：t，1ih10iihhtTsum?结束NoYes差分方程的解hik+1是否逼近原微分方程的解Hik+1?必须从差分方程的收敛性和稳定性两个方面回答此问题。显式差分格式收敛和稳定的条件2102)(Txt只有收敛和稳定的差分格式，才具有实用价值。因此，合理选取△x,△t是很重要的。收敛性：如果在△x,△t取得充分小时，差分方程的解和微分方程的解析解很接近，便说这种差分格式是收敛的。稳定性：差分计算时，每一步都有舍入误差，随着计算时间或计算次数的增加，累积误差逐渐减小，以至于不影响计算结果，那么这种差分格式便是稳定的。显式差分方程的收敛性和稳定性二、一维隐式有限差分格式tHxHT*22(i,k+1)(i-1,k+1)(i+1,k+1)(i,k)211111)(2xhhhTkikikithhkiki1*向后差分=kikikikikihhhhhxtT1111112*2)(整理得：截断误差为：O([Δx]2)+O(Δt)2*)(xtT定义kikikikihhhh11111)21(左端包含了三个未知数，不能直接解出hik+1，所以称为隐式差分格式。隐式差分方程的求解由于线性方程组的系数矩阵只在三条对角线上有值，其余均为零，所以称为三对角线方阵，可用追赶法求解。由于隐式差分格式不能直接解出hik+1值，所以应该对所有内结点（i=1,i=2,…,i=l-1)都按差分方程式列出相应的方程。并把边界条件代入形成由（l-1)个方程组成的线性方程组，联立求解。隐式差分格式是无条件收敛和稳定的。1122101112121112121212knknknkkknknkkhhhhhhhhhh开始输入初始参数t=0对i=0,1,…,n循环执行xiHhi00t=t+Δt对i=1,…,n-1循环线性方程组调用追赶法程序求解三对角线方程组输出结果：t，1ih10iihhtTsum?结束NoYeskikikikihhhh11111)21(三、一维中心式有限差分格式这种差分格式是在tk与tk+1时刻之间，取中间过渡时刻tk+△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([Δx]2)+O([Δt]2)21*kitHkikikixHxHTxHT2212221222121中心差分thhxhhhxhhhTkikikikikikikiki