计算土力学
主讲教师：张爱军
§ 2.7 岩土工程基本分析方法
? 分类
? 总应力分析法（ p＝ 0）
? 初期变形分析
使用不排水指标，计算出的结果是加荷瞬
间或短期的应力与变形。
?终期变形分析
采用排水指标，计算出的结果是最终或长
期的应力和变形，即：孔压消散完毕，主
固结完成后最终的应力和变形。
? 有效应力分析法（ p≠0）
§ 2.7.1 总应力分析
? 将土体视为固体，不考虑土中的渗流，不区分
单元中分别由土骨架和孔隙水压力传递和承受
的应力，也就是不考虑有效应力与孔隙水压力，
只考虑土体总应力。其分析基本理论与固体力
学中无任何区别
? 适用于不考虑渗流固结的情况，例如：饱和粗
粒土地基、强透水土料组成的土坝路堤的应力
和变形，以及饱和软粘土地基短期变形和稳定
问题。
? 其基本方程：
?平衡方程
?几何方程
?物理方程
?孔隙水渗流连续方程
?孔隙水平衡方程
我们以静力问题为例说明总应力分析
的基本方程。同时讲授基本方程的矩
阵表示方法。
? 平衡方程
其矩阵表示为：
0yxx zxx y z??????? ? ?? ? ?
0y x y y zx y z? ? ?? ? ??? ? ???
zyzx z
x y z g
?? ? ??? ???
? ? ??
? ? ? ? ? ?TBR? ??
? ? ? ? ? ?TBR? ?
0 0 0
0
0 0 0 0
0 0 0
00
00
00
0
0
0
x
y
z
xy
yz
zx
T
x
y
z
xy
yz
zx
x y z
y x z
g
z y z
x
y
z
yx
zy
zx
?
?
?
?
?
?
?
?
?
?
?
?
?
????
? ? ?
????
? ? ?
????
??
????
? ? ? ? ? ? ?
??
? ? ? ???
? ? ?
? ? ? ???
?
??
????
? ? ?
????
? ? ?
????
??
???
??
?
??
?
??
?? ?
?
??
?
??
??
?
?? ?
??
??
????
??
??
??
?? ??
??
??
??
??
????
0
0
g?
?
??
??
??
??
? ? ? ?
??
? ? ?
? ? ? ?
?
??
??
??
??
??
[B]T
[σ]
[σ]
? 物理方程
? 几何方程
? 总体控制方程就是：
将,与 代入：
得到：
? ? ? ?? ?D???
{ } [ ] { }
{ } {,,,,,,}
{ } {,,}
T
x y z x y y z x z
T
Bu
u u v w
?
? ? ? ? ? ? ?
??
其中：,为应变矩阵等于
：为位移矩阵等于
? ? ? ?? ?D??? ? ? ? ? ? ?TBR? ??? ? ? ?? ?Bu? ??
? ? ? ?? ?? ? ? ?TB D B u R?
? 这就是总应力法静力问题控制方程
这个方程实际上是 Biot动力固结方程中的惯性力
为 0，不考虑孔压的特殊形式。
? ? ? ?? ?? ? ? ?TB D B u R?
? ?? ? ? ?K u R?
? ? ? ? ? ? ? ? ? ?
? ?
? ?
? ? ? ? ? ? ? ?
T
K K B D B
B
D
S S D B
?
?
?
其中,－为刚度矩阵,＝
为几何矩阵
为物理矩阵
为应力矩阵,＝
进一步表示为：
? 对于各向异性、横观各向同性体其实就
是 [D]阵不一样而已，没有什么复杂的，
总体控制方程的矩阵形式是一致的。
? 但是对于 {R}阵，对于有面力的情况和有
初始应变或应力的情况，还有不同的形
式，以后的章节中再讲，但是整体控制
方程的矩阵形式是一样的。
? 对于动力情况以上控制方程是不适合的，
因为存在惯性力项，另外物理矩阵中还
应该包含粘性的成分，即：
? ? ? ? ? ?(,)f? ? ??
速度项
§ 2.7.2 有效应力分析
? 有效应力法要严格区分土骨架变形和孔
隙水压力，更能够真实地反映土体的自
身特性。
? 适用于：黏土地基、坝体和路基随时间
固结沉降分析，特别适合于求解随时间
变化的过程分析，反映施工过程。
? 其控制方程就是 Biot静力固结方程
?平衡方程
由于平衡方程只涉及总应力因此与总应
力法是一致的，即：
?几何方程，物理方程也是一致的，即：
?有效应力原理：表示成为矩阵形式：
? ? ? ? ? ?TBR? ??
? ? ? ?? ?D???? ? ? ?? ?,Bu? ??
? ? ? ? ? ?
? ? ? ?
'
1 1 1 0 0 0 T
Mp
M
?? ??
其中,＝
将物理方程、几何方程、有效应力原理
代入平衡方程中，得到：
?孔隙流平衡方程，若略去惯性力就为
Darcy定理，即：
? ? ? ? ? ? ? ? ? ? ? ? ? ?
? ?
TTB D B u B M p R
up
??
其中：未知量为,
wx
w
wy
w
wz
w
kp
v
gx
kp
v
gy
kp
v
gz
?
?
?
?
??
?
?
??
?
?
??
?
变成矩阵形式为：
?孔隙流体连续方程为：
{ } [ ] [ ] { }
{ } {,,}
00
1
[ ] 0 0
00
[ ] { }
T
w
T
w wx wy wz
v k B M p
v v v v
k
kk
g
k
BM
?
??
?
??
??
?
??
????
其中：
,与上含义相同 3× 3× 3× 6× 6
× 1＝ 3× 1
()wyw x w zvvv u v wx y z t x y z??? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ?
{ } [ ] { } ( { } { } )TTwM B v Mt ??? ?
x y z? ? ???
将 Darcy定理应用于连续方程中得到：
?这样就得到静力有效应力控制方程：
其中：未知量有,{u},p
还包含了时间项 t
{ } [ ] [ ] [ ] { } { } [ ] { } 0T T TM B k B M p M B ut????
? ? ? ? ? ? ? ? ? ? ? ? ? ?
{ } [ ] [ ] [ ] { } { } [ ] { } 0
TT
T T T
B D B u B M p R
M B k B M p M B u
t
??
?
??
?
§ 2.7.3 初始条件和边界条件
? 有了总体控制方程是无法求解的，还必须有初
始条件和边界条件
? 对于总应力法其未知量为三个方向的位移，因
此只涉及位移或应力初始条件和边界条件。对
于有效应力分析中还存在孔隙水压力未知量，
因此还增加渗透边界条件。
? 我们以有效应力分析为例说明边界条件和初始
条件问题，总应力分析中只要去掉渗透边界条
件即可。
? 边界条件
?位移边界条件
位移边界条件分为二种，即：约束边
界和已知位移边界。其中约束边界就
是在某个边界上某个自由度上位移受
到约束，在这个自由度上的位移为 0。
如图所示：
u=0 w=0v=0
——也可以是以上情况的组合 。x
y
z
例子：
或者：
另外一种是已知位移边界，即：
?应力边界
就是指在某个边界上存在面力。
?已知水压力边界，也就是在某个边界
上孔压或水头已知。
?透水边界：
?水头已知边界：
? ? ? ?uu?
0p?
,p p H H??或
?已知流量边界，也就是在某个边界上
渗透流量是已知的。
?不透水边界
?输入流量边界
设：边界上沿外法线方向的已知流速
为 vn，则流速边界条件可表示为：
L,m,n为边界面外法线方向与 xyz坐
标轴的夹角的余弦（方向余弦）。
x y z nlv mv nv v? ? ?
那么：根据达西定理：
对于不透水边界为：
对于已知流量边界为：
h h h
n
w w w
k k kp p pl m n v
g x g y g z? ? ?
? ? ?? ? ? ?
? ? ?
0h h h
w w w
k k kp p pl m n
g x g y g z? ? ?
? ? ?? ? ? ?
? ? ?
h h h
n
w w w
k k kp p pl m n v
g x g y g z? ? ?
? ? ?? ? ? ?
? ? ?
已知流速
也就是已知
流量
? 初始条件
?初始位移条件
?初始应力条件
?初始孔隙水压力或水头条件
?初始流量条件
这些条件是进行迭代计算的基础。
§ 2.7.4 有效应力法与总应力法
的关系
? 从控制方程可以看出：总应力法实际上是有效
应力法当 p＝ 0时的特殊形式，编制一个有效应
力程序是可以用作总应力计算的。
? 有效应力法是土力学区别于固体力学的主要方
面，作为一个土力学工作者，必须掌握。这也
是我们可以从事固体力学有些问题研究，而固
体力学人从事不了我们工作的主要原因。
? 总应力法作用也是很大的，也是目前工程上用
的较多的，也应该重视。
? 总应力法应用问题
?谢康和书中,TT”法，,Total stress
analysis with Terzaghi’s consolidation
theory” 。就是用总应力法采用排水指
标计算土体的最终变形和应力，然后
用 Terzaghi一维或二、三维固结理论
计算变形从开始到最终随时间变化的
过程。这个方法经过证明对于一维固
结是精确的，对于二三维固结是近似
的，但是也可以满足一定范围内工程
计算的需要。
?总应力法采用不排水指标可以计算土
体初期和瞬间变形和应力。
?总应力法采用排水指标可以计算土体
最终变形和应力。也就是孔隙水压力
消散完毕，主固结完成后最终的土体
变形和应力。
“UU,CU,CD； q,Cq,s,?
三轴试验
UU－不固结不排水试验
CU－固结不排水试验
CD－固结排水试验
直接剪切试验
q －快剪试验
Cq－固结快剪试验
s －慢剪试验
施工期或透水
小时
有先期固结应力
而不排水时，如：
地基深处的土体
运行期
?思考：在总应力分析中对于渗透水压
力该不该加，若该加那么如何施加？
? 至此 Biot固结理论的控制方程和初始、
边界条件均有了，并从数学理论上看，
未知数与方程的数量相同，方程是可解
的。
? 以下的工作有三个：
?寻求方程的求解方法
?寻求方程中参数的确定方法
?实现求解的手段分析，以及计算结果
的应用
第三章 有限差分法
? 基本思想：
用差分网格离散求解域，用差分公式将
控制方程（常微分方程或偏微分方程）
转化为差分方程（线性方程组），结合
边界条件和初始条件，求解线性方程组。
? 特点：
直观方便，编程容易。
? 发展历史
? 从 20世纪 40年代创立，发展至今
? 40年代后期，在土工渗流和固结问题分析中
成功应用
? 50年代及 60年代初，在弹性地基粱与板桩计
算中也应用有限差分法。
? 60年代后期，由于有限元法的发展，应用
较少。
? 目前由于 任意网格有限差分法 的出现，以及
与有限元完美结合使得这种方法成为一种普
遍的方法。
? 在求解 Biot固结方程中，空间域上的离散是
用有限元法，在时间域上的离散用差分法。
§ 3.1 差分公式
? 将求解域划分为差分网格（以二维问题
为例说明）
? 各点计算公式为：
边界节点
内部节点
外部虚拟节点
x
y
h
l
(xi,yi)
(x0,y0)
00,,iix x ih y y jl i j? ? ? ? 为 整 数
? 几个概念
? 区域：用 V表示；边界：用 S表示
? 步长,h,l分别为 x,y方向的步长
? 结点,两条平行线的交点，即（ xi,yi）
? 差分网格,x,y方向的平行线为差分网格，
若,h=l称正方形网格,h与 l不相等成为矩
形网格。
? 内部结点、边界结点，外部虚拟结点
处于 V中，处于 S上,处于 V外
? 将函数 f(x,y)在 结点 (xi,yi)上的值简记为 fi,j
? 导数的差分公式
?依据与函数的 Taylor展开式：
在点 (xi,yj)附近，二元函数 f(x,y)沿 x方
向的展开式为：
由上式：取 也就是
前后走一步，得：
23
23
,23
11(,) ( ) ( ) ( )
2 ! 3 !i j i i i
f f ff x y f x x x x x x
x x x
? ? ?? ? ? ? ? ? ?
? ? ?
,iix x h x x h? ? ? ?
2 2 3 3
1,,23
2 2 3 3
1,,23
2 ! 3!
2 ! 3!
i j i j
i j i j
f h f h f
f f h
x x x
f h f h f
f f h
x x x
?
?
? ? ?
? ? ? ?
? ? ?
? ? ?
? ? ? ?
? ? ?
当 h足够小时：可以忽略 3次幂以上各项，
得到简化计算式子：
联立求解得到：
这就是导数在 x方向上的差分公式
22
1,,2
22
1,,2
2!
2!
i j i j
i j i j
f h f
f f h
xx
f h f
f f h
xx
?
?
??
? ? ?
??
??
? ? ?
??
1,1,
2
1,,1,
22
2
2
i j i j
i j i j i j
fff
xh
f f ff
xh
??
??
??
?
?
???
?
?
同时，对于 y方向得到差分公式为：
还可以导出混合导数差分公式：
,1,1
2
,1,,1
22
2
2
i j i j
i j i j i j
fff
yl
f f ff
yl
??
??
??
?
?
???
?
?
2
,1,1
1,1 1,1 1,1 1,1
1,1 1,1 1,1 1,1
( ) ( )
2
[ ( ) ( ) ] /( 2 )
22
4
i j i j
i j i j i j i j
i j i j i j i j
ffff
x y x y x l
f f f f
h
ll
f f f f
hl
??
? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ?
?? ? ? ?
??
? ? ? ? ?
??
??
???
?
－－这就是中心差分公式
? 还有向前、向后差分公式，用得较少，
自己看一下。
? 对于不规则边界，可以采用以下差分公
式：
V
S
(i,j)
(i-1,j) (i+1,j)
(i,j-1)
(i,j+1)
(i,j+η)
(i+ξ,j)
h
ηl
ξh
l
,,1
2
,,1,22
( 1 )
1 2 2 2
[]
( 1 ) 1
i j i j
i j i j i j
fff
yl
f
f f f
yl
?
?
?
? ? ? ?
??
??
??
?
??
?
? ? ?
? ? ?
,1,
2
,1,,22
( 1 )
1 2 2 2
[]
( 1 ) 1
i j i j
i j i j i j
fff
xh
f
f f f
xh
?
?
?
? ? ? ?
??
??
??
?
??
?
? ? ?
? ? ?
Any questions？
主讲教师：张爱军
§ 2.7 岩土工程基本分析方法
? 分类
? 总应力分析法（ p＝ 0）
? 初期变形分析
使用不排水指标，计算出的结果是加荷瞬
间或短期的应力与变形。
?终期变形分析
采用排水指标，计算出的结果是最终或长
期的应力和变形，即：孔压消散完毕，主
固结完成后最终的应力和变形。
? 有效应力分析法（ p≠0）
§ 2.7.1 总应力分析
? 将土体视为固体，不考虑土中的渗流，不区分
单元中分别由土骨架和孔隙水压力传递和承受
的应力，也就是不考虑有效应力与孔隙水压力，
只考虑土体总应力。其分析基本理论与固体力
学中无任何区别
? 适用于不考虑渗流固结的情况，例如：饱和粗
粒土地基、强透水土料组成的土坝路堤的应力
和变形，以及饱和软粘土地基短期变形和稳定
问题。
? 其基本方程：
?平衡方程
?几何方程
?物理方程
?孔隙水渗流连续方程
?孔隙水平衡方程
我们以静力问题为例说明总应力分析
的基本方程。同时讲授基本方程的矩
阵表示方法。
? 平衡方程
其矩阵表示为：
0yxx zxx y z??????? ? ?? ? ?
0y x y y zx y z? ? ?? ? ??? ? ???
zyzx z
x y z g
?? ? ??? ???
? ? ??
? ? ? ? ? ?TBR? ??
? ? ? ? ? ?TBR? ?
0 0 0
0
0 0 0 0
0 0 0
00
00
00
0
0
0
x
y
z
xy
yz
zx
T
x
y
z
xy
yz
zx
x y z
y x z
g
z y z
x
y
z
yx
zy
zx
?
?
?
?
?
?
?
?
?
?
?
?
?
????
? ? ?
????
? ? ?
????
??
????
? ? ? ? ? ? ?
??
? ? ? ???
? ? ?
? ? ? ???
?
??
????
? ? ?
????
? ? ?
????
??
???
??
?
??
?
??
?? ?
?
??
?
??
??
?
?? ?
??
??
????
??
??
??
?? ??
??
??
??
??
????
0
0
g?
?
??
??
??
??
? ? ? ?
??
? ? ?
? ? ? ?
?
??
??
??
??
??
[B]T
[σ]
[σ]
? 物理方程
? 几何方程
? 总体控制方程就是：
将,与 代入：
得到：
? ? ? ?? ?D???
{ } [ ] { }
{ } {,,,,,,}
{ } {,,}
T
x y z x y y z x z
T
Bu
u u v w
?
? ? ? ? ? ? ?
??
其中：,为应变矩阵等于
：为位移矩阵等于
? ? ? ?? ?D??? ? ? ? ? ? ?TBR? ??? ? ? ?? ?Bu? ??
? ? ? ?? ?? ? ? ?TB D B u R?
? 这就是总应力法静力问题控制方程
这个方程实际上是 Biot动力固结方程中的惯性力
为 0，不考虑孔压的特殊形式。
? ? ? ?? ?? ? ? ?TB D B u R?
? ?? ? ? ?K u R?
? ? ? ? ? ? ? ? ? ?
? ?
? ?
? ? ? ? ? ? ? ?
T
K K B D B
B
D
S S D B
?
?
?
其中,－为刚度矩阵,＝
为几何矩阵
为物理矩阵
为应力矩阵,＝
进一步表示为：
? 对于各向异性、横观各向同性体其实就
是 [D]阵不一样而已，没有什么复杂的，
总体控制方程的矩阵形式是一致的。
? 但是对于 {R}阵，对于有面力的情况和有
初始应变或应力的情况，还有不同的形
式，以后的章节中再讲，但是整体控制
方程的矩阵形式是一样的。
? 对于动力情况以上控制方程是不适合的，
因为存在惯性力项，另外物理矩阵中还
应该包含粘性的成分，即：
? ? ? ? ? ?(,)f? ? ??
速度项
§ 2.7.2 有效应力分析
? 有效应力法要严格区分土骨架变形和孔
隙水压力，更能够真实地反映土体的自
身特性。
? 适用于：黏土地基、坝体和路基随时间
固结沉降分析，特别适合于求解随时间
变化的过程分析，反映施工过程。
? 其控制方程就是 Biot静力固结方程
?平衡方程
由于平衡方程只涉及总应力因此与总应
力法是一致的，即：
?几何方程，物理方程也是一致的，即：
?有效应力原理：表示成为矩阵形式：
? ? ? ? ? ?TBR? ??
? ? ? ?? ?D???? ? ? ?? ?,Bu? ??
? ? ? ? ? ?
? ? ? ?
'
1 1 1 0 0 0 T
Mp
M
?? ??
其中,＝
将物理方程、几何方程、有效应力原理
代入平衡方程中，得到：
?孔隙流平衡方程，若略去惯性力就为
Darcy定理，即：
? ? ? ? ? ? ? ? ? ? ? ? ? ?
? ?
TTB D B u B M p R
up
??
其中：未知量为,
wx
w
wy
w
wz
w
kp
v
gx
kp
v
gy
kp
v
gz
?
?
?
?
??
?
?
??
?
?
??
?
变成矩阵形式为：
?孔隙流体连续方程为：
{ } [ ] [ ] { }
{ } {,,}
00
1
[ ] 0 0
00
[ ] { }
T
w
T
w wx wy wz
v k B M p
v v v v
k
kk
g
k
BM
?
??
?
??
??
?
??
????
其中：
,与上含义相同 3× 3× 3× 6× 6
× 1＝ 3× 1
()wyw x w zvvv u v wx y z t x y z??? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ?
{ } [ ] { } ( { } { } )TTwM B v Mt ??? ?
x y z? ? ???
将 Darcy定理应用于连续方程中得到：
?这样就得到静力有效应力控制方程：
其中：未知量有,{u},p
还包含了时间项 t
{ } [ ] [ ] [ ] { } { } [ ] { } 0T T TM B k B M p M B ut????
? ? ? ? ? ? ? ? ? ? ? ? ? ?
{ } [ ] [ ] [ ] { } { } [ ] { } 0
TT
T T T
B D B u B M p R
M B k B M p M B u
t
??
?
??
?
§ 2.7.3 初始条件和边界条件
? 有了总体控制方程是无法求解的，还必须有初
始条件和边界条件
? 对于总应力法其未知量为三个方向的位移，因
此只涉及位移或应力初始条件和边界条件。对
于有效应力分析中还存在孔隙水压力未知量，
因此还增加渗透边界条件。
? 我们以有效应力分析为例说明边界条件和初始
条件问题，总应力分析中只要去掉渗透边界条
件即可。
? 边界条件
?位移边界条件
位移边界条件分为二种，即：约束边
界和已知位移边界。其中约束边界就
是在某个边界上某个自由度上位移受
到约束，在这个自由度上的位移为 0。
如图所示：
u=0 w=0v=0
——也可以是以上情况的组合 。x
y
z
例子：
或者：
另外一种是已知位移边界，即：
?应力边界
就是指在某个边界上存在面力。
?已知水压力边界，也就是在某个边界
上孔压或水头已知。
?透水边界：
?水头已知边界：
? ? ? ?uu?
0p?
,p p H H??或
?已知流量边界，也就是在某个边界上
渗透流量是已知的。
?不透水边界
?输入流量边界
设：边界上沿外法线方向的已知流速
为 vn，则流速边界条件可表示为：
L,m,n为边界面外法线方向与 xyz坐
标轴的夹角的余弦（方向余弦）。
x y z nlv mv nv v? ? ?
那么：根据达西定理：
对于不透水边界为：
对于已知流量边界为：
h h h
n
w w w
k k kp p pl m n v
g x g y g z? ? ?
? ? ?? ? ? ?
? ? ?
0h h h
w w w
k k kp p pl m n
g x g y g z? ? ?
? ? ?? ? ? ?
? ? ?
h h h
n
w w w
k k kp p pl m n v
g x g y g z? ? ?
? ? ?? ? ? ?
? ? ?
已知流速
也就是已知
流量
? 初始条件
?初始位移条件
?初始应力条件
?初始孔隙水压力或水头条件
?初始流量条件
这些条件是进行迭代计算的基础。
§ 2.7.4 有效应力法与总应力法
的关系
? 从控制方程可以看出：总应力法实际上是有效
应力法当 p＝ 0时的特殊形式，编制一个有效应
力程序是可以用作总应力计算的。
? 有效应力法是土力学区别于固体力学的主要方
面，作为一个土力学工作者，必须掌握。这也
是我们可以从事固体力学有些问题研究，而固
体力学人从事不了我们工作的主要原因。
? 总应力法作用也是很大的，也是目前工程上用
的较多的，也应该重视。
? 总应力法应用问题
?谢康和书中,TT”法，,Total stress
analysis with Terzaghi’s consolidation
theory” 。就是用总应力法采用排水指
标计算土体的最终变形和应力，然后
用 Terzaghi一维或二、三维固结理论
计算变形从开始到最终随时间变化的
过程。这个方法经过证明对于一维固
结是精确的，对于二三维固结是近似
的，但是也可以满足一定范围内工程
计算的需要。
?总应力法采用不排水指标可以计算土
体初期和瞬间变形和应力。
?总应力法采用排水指标可以计算土体
最终变形和应力。也就是孔隙水压力
消散完毕，主固结完成后最终的土体
变形和应力。
“UU,CU,CD； q,Cq,s,?
三轴试验
UU－不固结不排水试验
CU－固结不排水试验
CD－固结排水试验
直接剪切试验
q －快剪试验
Cq－固结快剪试验
s －慢剪试验
施工期或透水
小时
有先期固结应力
而不排水时，如：
地基深处的土体
运行期
?思考：在总应力分析中对于渗透水压
力该不该加，若该加那么如何施加？
? 至此 Biot固结理论的控制方程和初始、
边界条件均有了，并从数学理论上看，
未知数与方程的数量相同，方程是可解
的。
? 以下的工作有三个：
?寻求方程的求解方法
?寻求方程中参数的确定方法
?实现求解的手段分析，以及计算结果
的应用
第三章 有限差分法
? 基本思想：
用差分网格离散求解域，用差分公式将
控制方程（常微分方程或偏微分方程）
转化为差分方程（线性方程组），结合
边界条件和初始条件，求解线性方程组。
? 特点：
直观方便，编程容易。
? 发展历史
? 从 20世纪 40年代创立，发展至今
? 40年代后期，在土工渗流和固结问题分析中
成功应用
? 50年代及 60年代初，在弹性地基粱与板桩计
算中也应用有限差分法。
? 60年代后期，由于有限元法的发展，应用
较少。
? 目前由于 任意网格有限差分法 的出现，以及
与有限元完美结合使得这种方法成为一种普
遍的方法。
? 在求解 Biot固结方程中，空间域上的离散是
用有限元法，在时间域上的离散用差分法。
§ 3.1 差分公式
? 将求解域划分为差分网格（以二维问题
为例说明）
? 各点计算公式为：
边界节点
内部节点
外部虚拟节点
x
y
h
l
(xi,yi)
(x0,y0)
00,,iix x ih y y jl i j? ? ? ? 为 整 数
? 几个概念
? 区域：用 V表示；边界：用 S表示
? 步长,h,l分别为 x,y方向的步长
? 结点,两条平行线的交点，即（ xi,yi）
? 差分网格,x,y方向的平行线为差分网格，
若,h=l称正方形网格,h与 l不相等成为矩
形网格。
? 内部结点、边界结点，外部虚拟结点
处于 V中，处于 S上,处于 V外
? 将函数 f(x,y)在 结点 (xi,yi)上的值简记为 fi,j
? 导数的差分公式
?依据与函数的 Taylor展开式：
在点 (xi,yj)附近，二元函数 f(x,y)沿 x方
向的展开式为：
由上式：取 也就是
前后走一步，得：
23
23
,23
11(,) ( ) ( ) ( )
2 ! 3 !i j i i i
f f ff x y f x x x x x x
x x x
? ? ?? ? ? ? ? ? ?
? ? ?
,iix x h x x h? ? ? ?
2 2 3 3
1,,23
2 2 3 3
1,,23
2 ! 3!
2 ! 3!
i j i j
i j i j
f h f h f
f f h
x x x
f h f h f
f f h
x x x
?
?
? ? ?
? ? ? ?
? ? ?
? ? ?
? ? ? ?
? ? ?
当 h足够小时：可以忽略 3次幂以上各项，
得到简化计算式子：
联立求解得到：
这就是导数在 x方向上的差分公式
22
1,,2
22
1,,2
2!
2!
i j i j
i j i j
f h f
f f h
xx
f h f
f f h
xx
?
?
??
? ? ?
??
??
? ? ?
??
1,1,
2
1,,1,
22
2
2
i j i j
i j i j i j
fff
xh
f f ff
xh
??
??
??
?
?
???
?
?
同时，对于 y方向得到差分公式为：
还可以导出混合导数差分公式：
,1,1
2
,1,,1
22
2
2
i j i j
i j i j i j
fff
yl
f f ff
yl
??
??
??
?
?
???
?
?
2
,1,1
1,1 1,1 1,1 1,1
1,1 1,1 1,1 1,1
( ) ( )
2
[ ( ) ( ) ] /( 2 )
22
4
i j i j
i j i j i j i j
i j i j i j i j
ffff
x y x y x l
f f f f
h
ll
f f f f
hl
??
? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ?
?? ? ? ?
??
? ? ? ? ?
??
??
???
?
－－这就是中心差分公式
? 还有向前、向后差分公式，用得较少，
自己看一下。
? 对于不规则边界，可以采用以下差分公
式：
V
S
(i,j)
(i-1,j) (i+1,j)
(i,j-1)
(i,j+1)
(i,j+η)
(i+ξ,j)
h
ηl
ξh
l
,,1
2
,,1,22
( 1 )
1 2 2 2
[]
( 1 ) 1
i j i j
i j i j i j
fff
yl
f
f f f
yl
?
?
?
? ? ? ?
??
??
??
?
??
?
? ? ?
? ? ?
,1,
2
,1,,22
( 1 )
1 2 2 2
[]
( 1 ) 1
i j i j
i j i j i j
fff
xh
f
f f f
xh
?
?
?
? ? ? ?
??
??
??
?
??
?
? ? ?
? ? ?
Any questions？
