偏微分方程的离散化方法4

x
x
x i
x
2、 一阶后差商
P P(x) P( x x) ， P Pi Pi1
x
x
x i
x
3、 一阶中心差商
P P( x x) P(x x) ， P Pi1 Pi1
x
2x
x i
2x
忽略截断误差 O(x) 忽略截断误差 O(x) 忽略截断误差 O(x 2 )
到第 n 时刻的 P 值，也要用到第 n+1 时刻的 P 值。这种差分格式是隐式差分格 式。在隐式差分格式中：在点（i，n+1），用到（i-1，n+1）、（i+1，n+1）、（i，n） 三点。在隐式差分格式中：只有一个方程，2~3 个未知数，但稳定，精度好，广 泛使用。
以上方程的一般形式： ci Pi1 ai Pi bi Pi1 di ，形成三对角矩阵。
2!
3!
4!
(*)
P(x) xP(x) O(x)
P(x x) P(x) x P(x) (x / 2)2 P(x) (x / 2)3 P(x)
2
2
2!
3!
P(x) x P(x) O(x)
2
2
1、 一阶前差商
P P(x x) P( x) ， P Pi1 Pi
偏微分方程的 离散化方法
一、离散化的概念
油藏是非均质的，岩石和流体性质伴随时间常常是发生变化的，建立的偏微 分方程一般是非线性的，求解偏微分方程的解析解比较困难，常用数值求解。 目前工程上应用的离散化方法有：有限差分法、有限元法、边界元法、变分 法等。 离散化的核心是把整体分成若干单元来处理，而每个小单元的形状是规则的， 并可以认为是均质的，从而把形状不规则的非均质的问题转化为形状规则的 均质的问题——非线性问题线性化。 计算过程中可以控制精度。要求的精度越高，则需要划分的单元就越多，计 算工作量相应就越大，反之，单元划分得少些，计算工作量就小，但精度变 差些。 微分方程离散化，主要在空间和时间两方面被离散化
x x0
x
P lim P(x) P(x x)
x x0
x
P
P(x x) P(x x)
lim
x x0
2x
P
前差商 后差商 中心差商
x
函数 P(x+Δx)利用 Talor 公式逼近导数
P(x x) P(x) xP(x) x 2 P(x) x3 P(x) x 4 P (4) (x)

Pi n
t
P n1 i
(1 2 )Pin


(
Pn i1

Pn i1
)
，


t x
2
，截断误差：O(t
x2 )
从方程可以看出：如果已知第 n（本步时间）的值 Pin ，就可以求得第 n+1
时刻（下步时间）的值
P n1 i
。因此如初始条件，即
n=0
时各网格的
P
值已给定，
12
上式两端同除 x 2 ,整理得:
P '' (x) P(x x) 2P(x) P(x x) O(x 2 ) x 2
忽略二阶截断误差 O(x 2 )
2P x 2

P(x

x)

2P(x) x 2

P(x

x)
，
2P x 2

Pi1
2Pi x 2
y
x
无效网格 有效网格 点中心网格 块中心网格
z
x y
局部网格加密
模拟区网格图（井位、边界、断层）
五点法注水开发5年后XW3层含水饱和度分布图
五点法注水开发20年后XW3层含水饱和度分布图
z
r
混合网格
二、有限差分法----导数的差商逼近
P lim P(x x) P(x)

Pi
n 1,
j
)

n1/ 2
P P i1, j
n1/ 2 i 1, j
P n1 i, j

Pn i, j
2x 2
2x
t
四、边界条件的处理
（一）、内边界条件处理
定产条件：即井以一定产量 q 生产。如在网格（i，j）上有一口井，产量 q，
则可在渗流方程左边加上产量相，生产井 q 为负，注水井 q 为正。
（2）隐式差分：在点（i，j，n+1）的差分方程（图示）
P n1 i1, j

2
P n1 i, j

P n1 i1, j

P n1 i, j1

2
Pi
n1 ,j

P n1 i, j1

P n1 i, j

Pn i, j
x 2
y 2
t
若取正方形网格：则： x y
，
u

u(x, y, t) u(x x2 , y, t)
x x x1
x1
x x x1
x2
2
2

k u
k x x1 2
u(x

x1,
y, t) x1

u( x,
y,
t)

k x x2 2
u(x, y, t) u(x x2 , y, t) x2

2
Pi
n ,j
x 2
ห้องสมุดไป่ตู้
Pn i1, j

Pn i, j1

2
Pn i, j
y 2

Pn i, j1

P n1 i, j

Pn i, j
t
P n1 i, j

Pn i, j


(
Pi
n 1,
j

2Pi
n ,j

Pi
n 1,
j
)


(
Pn i, j
1

2
Pi
n ,j

Pn i, j
1
)


t
x 2
，

t
y 2
，截断误差： O(t
x2

y 2 )
该线性代数方程组在节点（i，j）列方程式，用到（i，j），（i+1，j），（i-1， j），（i，j+1），（i，j-1）五个点。 —显式：只有一个方程，1 个未知数，简单，精度较差，时间步长受到严格限制， 基本不用。
就可以依次求得以后各时间的 P 值。这种差分格式是显式差分格式。在显式差分
格式中：只有一个未知数 Pin1 ，由一个方程就可以求出。简单，精度较差，时间
步长受到严格限制，基本不用。
（2）隐式差分：利用 P（x，t）关于 t 的一阶向后差商和关于 x 的二阶差商， 在点（i，n+1）的差分方程：
P n1 i 1
P P( x x / 2) P( x x / 2) ，P Pi1/2 Pi1/2 忽略截断误差 O((x / 2)2 )
x
x
x i
x
1、 二阶差商
将方程(*)正负相加,可得: P(x x) P(x x) 2P(x) x 2 P '' (x) x 4 P (4) (x) .........
P n1 i, j1

P n1 i1, j

(4

1

)
Pi
n 1 ,j

P n1 i1, j

Pn i, j1

1

Pn i, j
该线性代数方程组在节点（i，j）列方程式，也要用到（i，j），（i+1，j），
（i-1，j），（i，j+1），（i，j-1）五个点，但时刻不同。
—隐式：只有一个方程，5 个未知数，稳定，广泛使用。
其一般形式是：
ci
,
j
Pi
n1 , j 1

ai,
j
P n1 i 1, j

ei,
j
P n1 i, j

bi,
j
P n1 i 1

di,
j
P n1 i. j 1

fi, j （形成五对
角矩阵）
五对角矩阵形式
1 2 3 4 512345123451234512345
1PP
P
2PPP
P
3
PPP
P
4
PPP
P
5
PP
P
1P
PP
P
2
P
PPP
P
3
P
PPP
P
4
P
PPP
P
5
P
PP
P
1
P
PP
P
2
P
PPP
P
3
P
PPP
P
4
P
PPP
P
5
P
PP
P
1
P
PP
P
2
P
PPP
P
3
P
PPP
P
4
P
PPP
P
5
P
PP
P
1
P
PP
2
P
PPP
3
P
PPP
4
P
PPP
5
P
PP
3、Crank_Nicolson 差分格式
Crank_Nicolson 差分格式（简称 C_N 格式）是综合显式和隐式格式而构建， 将空间二阶差商取为 n 时刻与 n+1 时刻的算术平均值，则：
Pi1
(用节点位置)
i
1、 一种常用二阶差商处理方法
k u k u
x
k
x



x x x1 2 x
x x x2 2
， x
1 2
(x1

x2
)
u
u(x x1, y, t) u(x, y, t)

Pn i 1,
j
(2
2x 2
t
)Pi,nj

Pn i1, j
截断误差：O(x2 t)
这种差分格式求解精度高，工作量与隐式差不多，在油藏数值模拟中经常采 用的格式之一。
4、其它差分格式
时间中心显式差分：
Pn i1, j

2
Pn i, j
x 2

Pn i 1,
j

P n1 i, j
定压条件：即井以一定流动压力 Pwf 生产，这时的 q 未知，可由给定的井
2Pin1 x 2

P n1 i1

P n1 i

Pi n
t
(1
2
) Pi n 1


(
P n1 i1

Pi
n 1 1
)

Pi n
从方程可以看出：如果已知第 n（本步时间）的值 Pin ，为了求得第 n+1 时刻（下
步时间）的值 Pin1 ，必须解一个线性代数方程组。即：要想求出 Pin1 值，需用
（1）离散空间：把所研究的空间划分成某种类型的网格， 大的空间转化为若干小单元组成，网格之间动态连接，通 常采用矩形网格（正方体）。
（2）离散时间：把研究的时间域分成若干小的时间段， 在每个时间段内，对问题求解，时间段之间有机连接。步 长大小取决于所要解决的实际问题。
离散空间
P
t
离散时间
1、网格系统 它有x，y两个自变量，在平面上用平行线分割成许多网格， 如考虑时间，则。编号：x→i，y→j，t→n。为步长（对三 维z→k）。 节点：网格的交点叫网格节点。取一些与边界s接近的网格 节点，把他们连成折线Sh，Sh所围成的区域记为Dh，Dh 内的节点为内部节点、边界上的节点为边界节点。
2、 等距网格就是指建立差分网格时，所采用的步长都是 相等的，反之称为不等距网格。
3、网格类型 常规网格系统： （1）块中心网格：用网格小块的几何中心来表示小块的坐标 （2）点中心网格：用节点的坐标来表示小块的坐标
块中心网格和点中心网格的离散点数不同，但最终形成一样的差分方程，只 有在处理边界条件时各有方便之处，块中心网格比较容易处理定流量边界， 点中心网格比较容易处理定压边界。 非常规网格系统： （1）局部网格加密 （2）混合网格 （3）多边形网格
三对角矩阵形式
1 2 3 4 512345
1PP
2PPP
3 PPP
4
PPP
5
P PP
1
PPP
2
PPP
3
PPP
4
PPP
5
PP
2、椭圆型方程： 二维不稳定渗流方程
2 P x 2

2P y 2

P t
采用：等距网格差分 （1）显示差分：在点（i，j，n）的差分方程（图示）
Pn i1, j
x x
x
Δx
Δx1 Δx2
三、有限差分方程的建立
1、抛物型方程：一维不稳定渗流方程：
2 P x 2

P t
（1）显示差分：利用 P（x，t）关于 t 的一阶向前差商和关于 x 的二阶差商，在 点（i，n）的差分方程。
Pn i 1
2Pin x 2

Pn i 1

P n1 i

P n1 i, j
2t
Dufort_Frankel 差分格式：
Pn i1, j

(
Pi
n 1 ,j

Pi
n ! ,j
)

Pn i1, j
x 2

P n1 i, j

P n1 i, j
2t
Douglas_Jones 校正差分格式：
（1）预报差分格式：
n1/ 2 P 2P P i1, j
1
(
Pi
n 1,
j
2

2
Pi
n ,j
x 2

Pn i1, j

P n1 i1, j
2Pi,nj1 x 2

P n1 i1, j
)

P n1 i, j

t
Pn i, j
整理：
P n1 i1, j
(2
2x 2
t
)Pi
n 1 ,j

P n1 i1, j

n1/ 2 i, j
n1/ 2 i1, j
P P n i1, j
n i1, j

P P n1/ 2
n
i, j
i, j
x 2
2x
t / 2
（2）修正差分格式：
(P n1 i1, j

2
P n1 i, j

P n1 i1, j
)

(
Pi
n 1,
j

2
Pi
n ,j