搜索
基于有限差分法的微分方程离散化求解【摘要】目前偏微分方程数值求解的方法主要有两种，即有限差分法和有限元方法。本文论述了基于有限差分法的微分方程求解，离散化过程，并对结果进行了分析。【关键词】有限差分法离散化数值模拟1.前言有限差分法是计算机数值模拟最早采用比较成熟的方法，该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域，是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表述简单。有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内必改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。用有限差分法求解偏微分方程必须把连续问题进行离散化，为此首先要对求解区域进行离散化。构造离散网格系统的目的在于将表现为非均系统的大尺度用若干可以近似为均匀系统的尺度（如网格）表征。构造差分形式就是对参数在一定的离散点中心网格或块中心网格上离散。其中，离散网格可以是空间离散网格，也可以是时间离散网格（即离散时间步长）。平面网格的形式是多种多样的，如矩形网格、柱形网格、多边形网格等。空间离散网格是和边界条件相关联的，一般来说，对于第一类边界条件采用点中心网格较方便，第二、三类边界条件采用块中心网格比较合适。2.微分方程的离散化2.1一阶偏导数的差商逼近设有函数u(x,y,t)，对其自变量x的偏导数可以表示成函数的差商的极限形式（1）在⑴式中，当自变量增量充分小时，如果能用比较简单的函数差商来代替偏导数，即（2）这样就可以把偏微分方程用差分方程代替，从而把难于求解的偏微分方程化成代数方程组。利用taylor级数可以说用(2)形式的差商来逼近一阶偏听偏信导数，其误差对δx来说是一阶的。式(2)是用前差商来代替一阶偏导数即（3）同理，后差商也可以用来代替一阶偏导数，且其误差也为ο（δx）。类似地，还可以用中心差商（4）来代替偏导数。经证明，用中心差商代替一阶偏导数，其误差为ο（δx2）。这就意味着用中心差商来逼近偏导数，当δx非常小时，是一种比前两种差商逼近更为精确的方法。但是，在具体求解一个微分方程时，不仅要考虑差商对导数的逼近精度，更主要的是把差分方程作为一个整体考虑。因此，油藏数值模拟中也不总是用它来逼近一阶导数。2.2二阶偏导数的差分逼近在渗流微分方程中，除含有求知函数的一阶偏导数外，还含有未知函数的二阶偏导数。二阶偏导数相当于对函数一阶偏导数。在x轴上取三点：x-δx、x、x+δx，如图1所示。令将在点x处展开为中心差商（5）再将分别在x-δx/2,x+δx/2处展成一阶中心差商，得（6）（7）将式（6）、（7）代入（5）得（8）式（8）即是以二阶差商来代替二阶偏导数的公式，这种形成的差商称为二阶中心差商。2.3抛物型方程的离散化方法在各类低渗透介质渗流问题中，大多属于不稳定渗流，描述这类渗流问题的数学方程不仅含有未知量对空间坐标的偏导数，而且含有各变量对时间坐标的偏导数。这类方程在数学上属抛物型方程。对这类问题求解时，不仅要求出各变量在空间的分布状态，而且还要求出这种分布在时间区域上的变化情况。对数值解法来说，就是要求出各变量在不同时刻在各空间点上的数值分布。因为这类方程所描述的是不稳定渗流过程，所以方程中的各个变量都是随时间而变的，而空间差分项本身并没有指明它应该在什么时间点上取值，因此对应于不同的时间取值方法就可以得到不同形式的差分方程。更主要的是，由于求解这类问题时得到的是一系列时间点上的参数分布，后一个时间点上的结果依赖于前一个时间点的结果，这就产生了求解误差在时间坐标上的传播与变化的问题，即所谓的解的稳定性问题。取最简单的一维方程进行讨论：（9）其初始条件为边界条件为式（9）左端是对空间坐标的偏导数，右端是对时间坐标的偏导数。为了对时间的偏导数离散化，取时间坐标上两个时刻的步长为δt。对方程⑼的右端以任意时刻tn的后差商进行近似，在任意点xi得到(10)整理(10)式，得(11)在(11)中，除了含有本节点的未知变量uin+1外，还含有相邻节点的未知变量和。因此，这样的方程在求解时，不能逐节点依次求解，而需要逐节点列成方程组联立求解，称为隐式格式。由于采用隐式格式时计算过程的每一时间步都要用联立方程组求解，所以它的工作量比显式格式明显增加。但由于它在解的稳定性方面具有显式格式不可比拟的优势，因此这种差分格式在油藏数值模拟中应用很广泛。一维隐式格式根据初始条件与边界条件最终可得到三对角线性方程。对于二维区域内抛物型方程隐式差分格式所形成的五对角方程。三维区域内形成一个七对角线性方程组。3.差分方程的相容性、收敛性、和稳定性（1）差分方程的相容性截断误差是微分方程在离散化而形成差分方程时所略去的余项，如果在i点上空间步长及时间步长都趋于零时，截断误差ri也能足够小而且趋于零，则此时的差分方程必然逼近原微分方程。这种逼近在数学上称为相容逼近，或者说差分方程和原微分方程是相容的。如果只考虑截断误差ri的整数阶数的话，则当ri=ο(hp)(p≥1)差分方程和原微分方程相容；当ri=ο(h0)差分方程和原微分方程不相容。这样就可以简单地根据误差的阶数来判断其相容性。（2）差分方程的稳定性为了讨论差分方程求解时的稳定性问题，首先对用关分方程时可能产生的误差作一些分析。以方程为例。设在一定的初始及边界条件下，方程在i点上的精确解为uin，而用差分方法求解方程，在计算过程中不引入任何误差的前提下，所得到的解为ui*n。由于差分方程只是微分方程的一种近似表达式，一般来，有uin≠ui*n称uin-ui*n为微分方程差分求解的整体截断误差。显然，整体截断误差的产生是由于用差分方程代替微分方程存在着截断误差而造成的。但是，在通常情况下，在进行计算时，无论采用电子计算机还是人工计算，所取的有数字都不可能是无限多。这样就不可避免地在计算过程中忽略掉有效数字位数以后的些数字，从而引入一定的误差，称为舍入误差。因此，解差分方程求出的解往往还不是ui*n。为了区别，将其记为uin。uin为差分方程的近似解。显然，舍入误差可表示为εin=ui*n-uin考虑了上述两种形式的误差后，得到微分方程差分求解的实际误差为：微分方程差分求解理论中的收敛性与稳定性是直接与上述两种形式的误差相联系的。收敛性要求在时间步长和空间步长趋于零时，整体截断误差也趋于零。而稳定性则要求在计算步数不断增多时，无论在计算开始或计算过程中所引入的任何误差都不会随计算过程而不断增大。在数学上，一个差分格式是稳定的，是指某计算阶段所产生的任何误差不在以后的进一步计算中被扩大。因为一旦这种误差能被不断扩大而成为越来越大的误差，必然使所得的结果变得上下摆动而毫无意义，这种差分格式就是不稳定的。如用公式表示，令εn为时刻tn所产生的误差，并假定在进一步的计算时不继续产生新的误差，则由于tn存在着误差εn，在tn+1时刻，误差变为εn+1，于是当时，差分格式是稳定；当时，差分格式是不稳定的。（3）差分方程的收敛性差分方程的收敛性，是指对于一个差分格式，如果在所考虑的区域内（包括空间坐标与时间坐标）的任一点上，当空间与时间步长趋于零时，差分方程的真解趋于原微分方程的真解，则称该差分格式是收敛的；否则称差分格式是不收敛的。显然，收敛性对于一个差分格式来说是必不可少的，一个不收敛的差分格式是没有任何意义的。对于差分方程的收敛性直接进行理论分析，通常要比稳定性分析困难得多。然而，由于收敛性、稳定性和相容性之间有比较密切的关系，可以根据相容性和稳定性的分析来推断其收敛性。差分方程分析中lax等价理论就阐明了这种关系：对一个相容逼近于原微分方程的差分方程来说，稳定性是收敛性的必要和充分条件，或者说收敛性和稳定性是等价的。这个理论有重要的理论价值，因为据此人们就可利用比较简单的相容性判别方法和稳定性分析方法来推断比较难以判别的收敛性，而不需要直接对它来进行论证。也就是说，只要证明了差分方程的稳定性和相容性，也就证明了它的收敛性。油藏数值模拟中常见到的线性抛物型方程，一般符合这个理论。对椭圆型方程差分格式，对物理上有意义的油藏模拟问题，其边值问题的解一般是光滑的，只要其差分方程满足相容性条件，而且也符合稳定性要求，则在一般情况下，当空间步长趋于零时，差分方程的解收敛于原微分方程的解。4结束语本文推导了有限差分方程的求解方法并对其相容性收敛性、和稳定性进行了分析，可以看出该方法大大减轻计算的压力，为了保证整个差分方程的隐式程度，各个求解变量的系数和导数也要参加迭代。如果导数项处理得不好，在迭代过程中就可能出现震荡现象，使计算过程无法进行。因此在实际的计算的过程中如何对系数进行更好的处理，即如何更好、有效的利用离散化方法都是我们在今后的研究中值得思考的问题。参考文献：［1］顾尔祚.流体力学有限差分法基础.上海.上海交通大学出版社，1988.6.［2］郭永存，卢德唐，马凌宵.低渗透油藏渗流的差分数值模拟[j].水动力学研究与进展，2004，19（3）：288-293.［3］骆祖江等.水、气二相渗流耦合模型全隐式联立求解[j].煤田地质与勘探，2001，10（6）.［4］donadw.peaceman.fundamentalsofnumericalreservoirsimulation.elsevierscientificpublishingcompany,1997.