偏微分方程的离散化方法研究

偏微分方程的离散化方法一、离散化的概念油藏是非均质的，岩石和流体性质伴随时间常常是发生变化的，建立的偏微分方程一般是非线性的，求解偏微分方程的解析解比较困难，常用数值求解。目前工程上应用的离散化方法有：有限差分法、有限元法、边界元法、变分法等。离散化的核心是把整体分成若干单元来处理，而每个小单元的形状是规则的，并可以认为是均质的，从而把形状不规则的非均质的问题转化为形状规则的均质的问题——非线性问题线性化。计算过程中可以控制精度。要求的精度越高，则需要划分的单元就越多，计算工作量相应就越大，反之，单元划分得少些，计算工作量就小，但精度变差些。微分方程离散化，主要在空间和时间两方面被离散化（1）离散空间：把所研究的空间划分成某种类型的网格，大的空间转化为若干小单元组成，网格之间动态连接，通常采用矩形网格（正方体）。（2）离散时间：把研究的时间域分成若干小的时间段，在每个时间段内，对问题求解，时间段之间有机连接。步长大小取决于所要解决的实际问题。离散空间tP离散时间1、网格系统它有x，y两个自变量，在平面上用平行线分割成许多网格，如考虑时间，则。编号：x→i，y→j，t→n。为步长（对三维z→k）。节点：网格的交点叫网格节点。取一些与边界s接近的网格节点，把他们连成折线Sh，Sh所围成的区域记为Dh，Dh内的节点为内部节点、边界上的节点为边界节点。2、等距网格就是指建立差分网格时，所采用的步长都是相等的，反之称为不等距网格。3、网格类型常规网格系统：（1）块中心网格：用网格小块的几何中心来表示小块的坐标（2）点中心网格：用节点的坐标来表示小块的坐标块中心网格和点中心网格的离散点数不同，但最终形成一样的差分方程，只有在处理边界条件时各有方便之处，块中心网格比较容易处理定流量边界，点中心网格比较容易处理定压边界。非常规网格系统：（1）局部网格加密（2）混合网格（3）多边形网格无效网格有效网格点中心网格块中心网格xyyzx局部网格加密模拟区网格图（井位、边界、断层）五点法注水开发5年后XW3层含水饱和度分布图五点法注水开发20年后XW3层含水饱和度分布图rz混合网格二、有限差分法----导数的差商逼近xxPxxPxPx)()(lim0前差商xxxPxPxPx)()(lim0后差商xxxPxxPxPx2)()(lim0中心差商Px函数P(x+Δx)利用Talor公式逼近导数)()()()(!4)(!3)(!2)()()()4(432xOxPxxPxPxxPxxPxxPxxPxxP(*))2()(2)()(!3)2/()(!2)2/()(2)()2(32xOxPxxPxPxxPxxPxxPxxP1、一阶前差商xxPxxPxP)()(，xPPxPiii1忽略截断误差)(xO2、一阶后差商xxxPxPxP)()(，xPPxPiii1忽略截断误差)(xO3、一阶中心差商xxxPxxPxP2)()(，xPPxPiii211忽略截断误差)(2xOxxxPxxPxP)2/()2/(，xPPxPiii2/12/1忽略截断误差))2/((2xO1、二阶差商将方程(*)正负相加,可得:.........)(12)()(2)()()4(4''2xPxxPxxPxxPxxP上式两端同除2x,整理得:)()()(2)()(22''xOxxxPxPxxPxP忽略二阶截断误差)(2xO222)()(2)(xxxPxPxxPxP，211222xPPPxPiiii(用节点位置)1、一种常用二阶差商处理方法xxukxukxkxxxxx2221，)(2121xxx112),,(),,(1xtyxutyxxuxuxx，222),,(),,(1xtyxxutyxuxuxxxxtyxxutyxukxtyxutyxxukxukxxxxx222112),,(),,(),,(),,(21Δx1Δx2Δx三、有限差分方程的建立1、抛物型方程：一维不稳定渗流方程：tPxP22（1）显示差分：利用P（x，t）关于t的一阶向前差商和关于x的二阶差商，在点（i，n）的差分方程。tPPxPPPninininini12112)()21(111ninininiPPPP，2xt，截断误差：)(2xtO从方程可以看出：如果已知第n（本步时间）的值niP，就可以求得第n+1时刻（下步时间）的值1niP。因此如初始条件，即n=0时各网格的P值已给定，就可以依次求得以后各时间的P值。这种差分格式是显式差分格式。在显式差分格式中：只有一个未知数1niP，由一个方程就可以求出。简单，精度较差，时间步长受到严格限制，基本不用。（2）隐式差分：利用P（x，t）关于t的一阶向后差商和关于x的二阶差商，在点（i，n+1）的差分方程：tPPxPPPninininini12111112ninininiPPPP)()21(11111从方程可以看出：如果已知第n（本步时间）的值niP，为了求得第n+1时刻（下步时间）的值1niP，必须解一个线性代数方程组。即：要想求出1niP值，需用到第n时刻的P值，也要用到第n+1时刻的P值。这种差分格式是隐式差分格式。在隐式差分格式中：在点（i，n+1），用到（i-1，n+1）、（i+1，n+1）、（i，n）三点。在隐式差分格式中：只有一个方程，2~3个未知数，但稳定，精度好，广泛使用。以上方程的一般形式：iiiiiiidPbPaPc11，形成三对角矩阵。三对角矩阵形式12345123451PP2PPP3PPP4PPP5PPP1PPP2PPP3PPP4PPP5PP2、椭圆型方程：二维不稳定渗流方程tPyPxP2222采用：等距网格差分（1）显示差分：在点（i，j，n）的差分方程（图示）tPPyPPPxPPPnjinjinjinjinjinjinjinji,1,21,,1,2,1,,122)2()2(1,,1,,1,,1,1,njinjinjinjinjinjinjinjiPPPPPPPP2xt，2yt，截断误差：)(22yxtO该线性代数方程组在节点（i，j）列方程式，用到（i，j），（i+1，j），（i-1，j），（i，j+1），（i，j-1）五个点。—显式：只有一个方程，1个未知数，简单，精度较差，时间步长受到严格限制，基本不用。（2）隐式差分：在点（i，j，n+1）的差分方程（图示）tPPyPPPxPPPnjinjinjinjinjinjinjinji,1,211,1,11,21,11,1,122若取正方形网格：则：yxnjinjinjinjinjinjiPPPPPP,1,1,11,1,111,1)14(该线性代数方程组在节点（i，j）列方程式，也要用到（i，j），（i+1，j），（i-1，j），（i，j+1），（i，j-1）五个点，但时刻不同。—隐式：只有一个方程，5个未知数，稳定，广泛使用。其一般形式是：jinjijinijinjijinjijinjijifPdPbPePaPc,11.,11,1,,1,1,11,,（形成五对角矩阵）五对角矩阵形式12345123451234512345123451PPP2PPPP3PPPP4PPPP5PPP1PPPP2PPPPP3PPPPP4PPPPP5PPPP1PPPP2PPPPP3PPPPP4PPPPP5PPPP1PPPP2PPPPP3PPPPP4PPPPP5PPPP1PPP2PPPP3PPPP4PPPP5PPP3、Crank_Nicolson差分格式Crank_Nicolson差分格式（简称C_N格式）是综合显式和隐式格式而构建，将空间二阶差商取为n时刻与n+1时刻的算术平均值，则：tPPxPPPxPPPnjinjinjinjinjinjinjinji,1,21,11,1,12,1,,1)22(21整理：njinjinjinjinjinjiPPtxPPPtxP,1,2,11,11,21,1)22()22(截断误差：)(2txO这种差分格式求解精度高，工作量与隐式差不多，在油藏数值模拟中经常采用的格式之一。4、其它差分格式时间中心显式差分：tPPxPPPnjinjinjinjinji221,1,2,1,,1Dufort_Frankel差分格式：tPPxPPPPnjinjinjinjinjinji2)(1,1,2,1!,1,,1Douglas_Jones校正差分格式：（1）预报差分格式：2/22,2/1,,1,122/1,12/1,2/1,1tPPxPPxPPPnjinjinjinjinjinjinji（2）修正差分格式：tPPxPPxPPPPPPnjinjinjinjinjinjinjinjinjinji,1,2/1,12/1,12,1,,11,11,1,122)2()2(四、边界条件的处理定产条件：即井以一定产量q生产。如在网格（i，j）上有一口井，产量q，则可在渗流方程左边加上产量相，生产井q为负，注水井q为正。定压条件：即井以一定流动压力Pwf生产，这时的q未知，可由给定的井底流动压力Pwf和井点所在网格节点的压力Pij计算：把网格内井近似看成稳态流动，符合平面径向流：)/ln()(2,,wewfjijirrPPKhq如考虑油层弹性和油井表皮效应的影响，则：)43ln()(2,,srrPPKhqwewfjiji（一）、内边界条件处理关于供给半径re,不同的学者,有不同的公式：1）、对于各向同性地层，步长：yx,，其等值供给半径：2214.0yxryxree2）、对于各向异性地层，步长：yx,，其等值供给半径：4/14/12/122/122/128.0yxxyyxxyeKKKKyKKxKKr（二）、外边界条件处理封闭边界：常取块中心网格并在边界网格外虚拟一排网格，并令其相邻两个网格压力相等。定压边界：常取点中心网格，由于边界点的压力一定，因此，只需求内部节点压力。