有限差分法的介绍及简单应用

工程中常用得数值解法有有限单元法和 差分法。 有限单元法 是以有限个单元的集合体来 代替连续体，属于物理上的近似。 差分法 是把弹性力学的基本方程和边 界条件（一般均为微分方程）近似地改用差 分方程（代数方程）来表示，把求解微分方 程的问题改换成为求解代数方程的问题，属 于数学上的近似。
我们在弹性体上，用相隔等 间距h而平行于坐标轴的两组平行 线织成正方形网格，Δx=Δy=h， 如图。 设f=f(x,y)为弹性体内的某一 个连续函数。该函数在平行于x 轴的一根网线上，如在３-０-１ 上,它只随x坐标的改变而变化。 在邻近结点０处，函数f可展为 泰勒级数如下：
(1 − 5)
差分公式(１-１)及(１-３)是以相隔2h的两 结点处的函数值来表示中间结点处的一阶导数 值，可称为中点导数公式。 以相邻三结点处的函数值来表示一个端点 处的一阶导数值，可称为端点导数公式。 应当指出：中点导数公式与端点导数公式 相比，精度较高。因为前者反映了结点两边的 函数变化，而后者却只反映了结点一边的函数 变化。因此，我们总是尽可能应用前者，而只 有在无法应用前者时才不得不应用后者。
至此，我们解决了怎样 计算边界上各结点
∂ϕ ∂ϕ ϕ, , ∂x ∂y
的值的问题。 至于边界外一行虚结点处 ϕ 的值，则可用边界上 结点处的 ∂ϕ 或
∂x
值来表示。例如，对于图1中的虚结点14，因为有
ϕ14 − ϕ 9 ∂ϕ = 2h ∂x 13
∂ϕ 值和边界内一行相应结点处 ∂y
ϕ 的
所以有
∂ϕ ϕ14 = ϕ 9 + 2h ∂x 13
当求出全部结点上的φ值以后，我们就可按应 力分量的差分公式（2－1）计算应力分量。 用差分法解弹性平面问题时，可按下列步骤进 行： （1） 在边界上任意选定一个结点作为基点A， 取
ϕA =
∂ϕ ∂ϕ =0 = ∂x A ∂y A
（5）按照公式（2－1）计算应力的分量。 说明： 1．以上是针对单连体导出的结果。对于多连 体，情况就不象这样简单。 2. 如果一部分边界是曲线的，或是不与坐标 轴正交，则边界附近将出现不规则的内结点。对 于这样的结点，差分方程（2－2）必须加以修正。
整理即得
20 0 −8(ϕ1 +ϕ2 +ϕ3 +ϕ4) + 2(ϕ5 +ϕ6 +ϕ7 +ϕ8) +(ϕ9 +ϕ10 +ϕ11 +ϕ12) = 0 ϕ (2−2)
对于弹性体边界以内的每一结点，都可以建立 这样一个差分方程。
二、联立求解这些线性代数方程，就 能求得各内结点处的值。 一般建立和求解差分方程，在数学上不 会遇到很大困难。但是，当对于边界内一行 的（距边界为h的）结点，建立的差分方程 还将涉及边界上各结点处的φ值，并包含边 界外一行的虚结点处的φ值。 为了求得边界上各结点处的φ值，须要 应用应力边界条件，即：
第二节
应力函数的差分解
当不计体力时，我们 已把弹性力学平面问题归 结为在给定边界条件下求 解双调和方程的问题。用 差分法解平面问题，就应 先将双调和方程变换为差 分方程，而后求解之。
一旦求得弹性体全部节点的φ值后，就可按 应力分量差分公式（对节点0）算得弹性体各节 点的应力。
(σ x )0
∂ 2ϕ = ∂y 2 ∂ 2ϕ = ∂x 2 1 = 2 [( ϕ 2 + ϕ 4 ) − 2ϕ 0 ] 0 h 1 = 2 [( ϕ 1 + ϕ 3 ) − 2ϕ 0 ] 0 h ( 2 − 1)
(1 − 1)
(1 − 2)
同理，在网线４-０-２上可得到差分公式
∂f ∂y f − f4 = 2 2h 0 f + f4 − 2 f0 = 2 2h 0 (1 − 3 ) (1 − 4 )
∂2 f ∂y 2
以上(１-１)~(１４)是基本差分公式,从 而可导出其它的差分 公式如下：
(d )
由式(d)及式(c)可见，设 即可根据面力分量及求得
ϕ A ,
∂ϕ ∂ϕ , , ∂x A ∂y A
为已知，
ϕ B ,
∂ϕ ∂ϕ , ∂y . ∂x B B
由前知，把应力函数加上一个线性函数，并不影 响应力。因此，可设想把应力函数加上a+bx+cy， 然后调整a，b，c三个数值，使得
由此பைடு நூலகம்：
d ∂ ϕ d s ∂y
= px = py
d ∂ ϕ − d s ∂x
关于边界上任一点处
∂ϕ ∂x
∂ϕ ∂y
∂ϕ ∂y = A
B B
的值，可将上式从A点到B点对s 积分得到：
∫
B
A
px d s
B
∂ϕ − = ∂x A
第一节
差分方程
从弹性力学的基本方程建立以来，这些方 程在各种问题的边界条件下如何求解，一直是 很多数学工作者和力学工作者研究的内容。即 弹性力学的经典解法存在一定的局限性，当弹 性体的边界条件和受载情况复杂一点，往往无 法求得偏微分方程的边值问题的解析解，许多 工程重要问题，不能够得出函数式的解答。因 此，弹性力学问题的各种数值解法便具有重要 的实际意义。 差分法是沿用已久的一种数值解法。随着 计算机的普及和相应的软件发展，此法成为解 弹性力学问题的一种有效的方法。
ϕA = 0
∂ϕ ∂ϕ = 0, ∂y = 0 ∂x A A
于是式(d)，式(c)即简化为：
∂ϕ ∂y ∂ϕ ∂x
B = ∫ px d s A B B = −∫ p y d s A B B B
( 2 − 3) ( 2 − 4) (2 − 5)
然后由面力的矩及面力之和算出边界上所有各 ∂ϕ ϕ 结点处φ 的值，以及所必需的一些 及 ∂值，即垂 直于边界方向的导数值。
∂x
∂y
（2）应用公式（2－6），将边界外一行虚结 点处的φ值用边界内的相应结点处的φ值来表示。 （3）对边界内的各结点建立差分方程（2－2）， 联立求解这些结点处的值。 （4）按照公式（2－6），算出边界外一行的各 虚结点处的φ值。
h2 ∂2 f ∂f f1 = f 0 + h + 2 ∂x 2 ∂x 0
0
0
(c )
(d )
联立(c),(d),解得差分公式：
f − f3 ∂f = 1 2h ∂x 0
∂2 f 2 ∂x f + f3 − 2 f0 = 1 h2 0
∂ 2ϕ ∂ 2ϕ − l1 ∂x∂y + l2 ∂x 2 = p y
由图2可见， l1=cos(N,x)=cosα=dy/ds, l2=cos(N,y)=sinα=-dx/ds, 于是，式（a）可改写为：
d y ∂ 2ϕ 2 d s ∂y d x ∂ 2ϕ + d s ∂x∂y = p x
将此式亦从A点到B点沿s进行积分，就得到边 界上任一点B处的φ值。为此利用分部积分法，得：
(ϕ )
B A
∂ϕ = x − ∂x A
B
∫
d ∂ϕ x A d s ∂x
B
∂ϕ d s + y ∂y − A
B
∫
B
A
y
d ∂ϕ d s ∂y
∂2 f 1 = 2 [( f 6 + f8 ) − ( f 5 + f 7 )] ∂x∂y 0 4h ∂4 f 1 4 = 4 [6 f 0 − 4( f1 + f 3 ) + ( f 9 + f11)] ∂x h 0 ∂4 f 1 2 2 = 4 [4 f 0 − 2( f1 + f 2 + f 3 + f 4 ) + ( f 5 + f 6 + f 7 + f8 )] ∂x ∂y h 0 ∂4 f 1 4 = 4 [6 f 0 − 4( f 2 + f 4 ) + ( f10 + f12 )] ∂y h 0 (1 − 6)
1 ∂2 f ∂f f = f 0 + ( x − x0 ) + 2 2! ∂x ∂x 0 1 ∂3 f 2 ( x − x0 ) + 3 3! ∂x 0 ( x − x0 )3 + ... 0
我们将只考虑离开结点０充分近的那些结 点，即（x-x0 ）充分小。于是可不计（x-x0 ）的 三次及更高次幂的各项，则上式简写为：
d s,
将式(b),(c)代入，整理得：
ϕ B = ϕ A + ( xB − xA )
B B ∂ϕ ∂ϕ + ( yB − y A ) + ∫ ( yB − y) px d s + ∫ ( x − xB ) py d s ∂y A ∂x A A A
l1σ x + l 2τ xy = p x l1τ xy + l 2σ
y
= py
∂ 2ϕ ∂ 2ϕ ∂ 2ϕ σ x = 2 , σ y = 2 , τ xy = − ∂y ∂x ∂x∂y
代入上式，即得：
∂ 2ϕ l1 2 ∂y ∂ 2ϕ − l2 ∂x∂y = p x
∫
A
py d s
∂ϕ ∂y ∂ϕ ∂x
∂ϕ = ∂y + B A ∂ϕ = − B ∂x A
∫ ∫
B
A B
px d s py d s
A
由高等数学可知，
d ϕ ∂ϕ d x ∂ϕ d y = + . . d s ∂x d s ∂y d s
d y ∂ 2ϕ d x ∂ 2ϕ − ∂x∂y − d s ∂x 2 = p y ds
d y ∂ 2ϕ d x ∂ 2ϕ 2 + ∂y d s ∂x∂y = p x ds d y ∂ 2ϕ d x ∂ 2ϕ − ∂x∂y − d s ∂x 2 = p y ds
第十章 有限差分法 弹性力学的经典解法存在一定的局限性，当 弹性体的边界条件和受载情况复杂一点，往往无 法求得偏微分方程的边值问题的解析解。因此， 各种数值解法便具有重要的实际意义。 常用数值解法有有限单元法和差分法。 第一节 第二节 差分公式的导出 应力函数的差分解
第三节 例 深梁的应力函数差分解
(σ )
y 0
(τ )
xy 0
∂ 2ϕ 1 = = − [( ϕ 5 + ϕ 7 ) − (ϕ 6 + ϕ 8 )] ∂x∂y 4h 2 0
可见，用差分法解平面 问题，共有两大任务： 一、建立差分方程 将（1-6）代入双调和方程
∂ 4ϕ ∂ 4ϕ ∂ 4ϕ +2 2 2 + =0 4 4 ∂x ∂x ∂y ∂y
1 ∂2 f ∂f f = f0 + ( x − x0 ) + 2 2! ∂x ∂x 0 ( x − x0 )2 0 (b)
在结点３，x=x0-h, 在结点1, x=x0+h，代入(b) 得： x=x 1, +h (b)
h2 ∂2 f ∂f f 3 = f 0 − h + 2 ∂x 2 ∂x 0
ϕ B = ∫ ( yB − y ) p x d s + ∫ ( x − xB ) p y d s
A A
从图易看出，式（2－3）右 边的积分式表示A与B之间的，x 方向的面力之和；式（2－4）右 边的积分式表示A与B之间的，y 方向的面力之和；式（2－5）右 边的积分式表示A与B之间的面 力对于B点的矩。