有限差分法及其应用

1、有限差分法及其应用1有限差分法简介有限差分法（FDM ）是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方程将解域 划分为差分网格，用有限个网络节点代替连续的求解域。有限差分法通过泰勒级数展开等方 法，把控制方程中的导数用网格节点上的函数值得差商代替进行离散，从而建立以网格节点 上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似值解法， 数学概念直观，表达简单，是发展较早且比较成熟的数值方法。2有限差分法的数学基础有限差分法的数学基础是用差分代替微分，用差商代替微商而用差商代替微商的意义是用函 数在某区域内的平均变化率来代替函数的真是变化率。而根据泰勒级数展开可以看出

2、，用差 商代替微商必然会带来阶段误差，相应的用差分方程代替微分方程也会带来误差，因此，在 应用有限差分法进行计算的时候，必须注意差分方程的形式，建立方法及由此产生的误差。3有限差分解题基本步骤有限差分法的主要解题步骤如下：1）建立微分方程根据问题的性质选择计算区域，建立微分方程式，写出初始条件和边界条件。2）构建差分格式首先对求解域进行离散化，确定计算节点，选择网格布局，差分形式和步长；然后以有 限差分代替无线微分，以差商代替微商，以差分方程代替微分方程及边界条件。3）求解差分方程差分方程通常是一组数量较多的线性代数方程，其求解方法主要包括两种:精确法和近似 法。其中精确法又称直接发，主要包括

3、矩阵法，高斯消元法及主元素消元法等；近似法 又称间接法，以迭代法为主，主要包括直接迭代法，间接迭代法以及超松弛迭代法。4）精度分析和检验对所得到的数值进行精度与收敛性分析和检验。4商用有限差分软件简介商用有限差分软件主要包括FLAC、UDEC/3DEC和PFC程序，其中，FLAC是一个基于显式 有限差分法的连续介质程序，主要用来进行土质、岩石和其他材料的三维结构受力特性模拟 和塑性流动分析；UDEC/3DEC是针对岩体不连续问题开发，用于模拟非连续介质在静，动态 载荷作用下的反应；PFC是利用显式差分算法和离散元理论开发的微、细观力学程序，它是 从介质的基本粒子结构的角度考虑介质的基本力学特性

4、，并认为给定介质在不同应力条件下 的基本特征主要取决于粒子之间接粗状态的变化，适用于研究粒状集合体的破裂和破裂发展 问题，以及颗粒的流动（大位移）问题。5有限差分法的应用1）光子晶体光纤色散的有限差分法研究光子晶体光纤的横截面结构和特征参数如图1所示。图1光子品体光纤横截面结构及特征参数示意图对结构如此复杂的光子晶体光纤，采用有限差分法加以研究。有限差分法是光波导分析中广 为采用的数值方法之一，它利用台劳展开用有限差分式代替波动方程中的微分式得到关于场 分量的有限差分方程。在纵向（z方向）均匀分布的光子晶体光纤中，电场的横向分量Ex和Ey 满足如下耦合方程组（卮 &一时&备J = 0（1）J

5、E. , ） 1 强 L E I 于+号铲+ W +（也,一 口旗,.+瑚;鼻志E）= 0（2）其中,Ex和Ey为电场的两个偏振分量,8 r为相对介电常数,k0为真空中的波数,B =koneff为传 播常数,neff为模式有效折射率。在光子晶体光纤分析中,通常求解基模（类比于传统光纤，也称 为HE11模)的两个偏振分量。文献中已经证明了横向结构对称分布的光子晶体光纤基模的线 偏振特性，因而(1),式中的电场偏振分量Ex和Ey可以认为不再耦合此时(1),式将变为两 个独立的分别关于两个偏振分量Ex和Ey的半矢量形式的波动方程3. 3 )史工=03)以捐8 0(4)(3),(4)式中保留左边的第二

6、项对于研究空气填充率较高的光子晶体光纤是必要的，因为这时 标量波动方程给出的解会产生较大的误差。在x,y方向采用相同的网格间距h,利用五点差分 格式可以建立基于半矢量波动方程(3),(4)的差分方程。J E夕I. +展土却+ %已中叩I +您E51 +(% + % )Ez + 篇 Q (p，0 ) 8 j Ez = 0(5)其中E为Ex或Ey,p,q为网格点编号。对于Ex,各项系数分别为=2 X L .)%- h 9 3 * 0)+ 豚(P h 弓)=_2& (p+ L g)% = j + . + 弟=& %对于Ey-各项系数分别为(5)式可化为特征值方程tI E)=矿 E=好蓦 E(6 )其

7、中A为系数矩阵,E为由各网格点构成的向量。计算结果为将基于半矢量波动方程的差分法得到的结果与文献中全矢量方法得到的结果比较，取与文 献中相同的求解参数:光纤参数为空气孔间距A=2.3 u m,空气孔直径d=1.0u m,整个求解区域 取为6AX6A,石英的折射率为1.45,波长为1.5 um,采用电壁(electric wall)边界条件,x,y方向 格点数均为240。基于半矢量波动方程的差分法得到结果为neff=1.42806，与文献中给出的结 果1.42868以及其他方法得到的结果较为接近。考虑到全矢量方法计算结果表明在求解某一 偏振态时另一偏振分量并不为零，但相差几个数量级，而半矢量波动

8、方程完全忽略了另一偏振 分量的影响，以及网格剖分的方法不同这一结果是较为精确的。2)抗滑桩全桩内力计算“m-k”法的有限差分法如图1所示，设抗滑桩全长为H,其中滑动而以上，即受荷段桩长为h1;滑动而以下， 即锚固段桩长为h2。设滑坡推力ET按梯形分布，地而和滑动而处的分布荷载集度分别为q0 和qa。桩前滑体的剩余抗滑力为E。桩的截而尺寸:长x宽为axb，桩的抗弯刚度为EI。为便 于公式的推导，设桩前滑体地基系数为K1 (z )=mZ,滑而以下锚固体地基系数为K2 (z ) = K常 数)。图1淖悭抗滑桩受力和变位Fig.l Forces acting Bit elastic stabilizi

9、ngpil* and its disphcemenc由于抗滑桩有刚性桩和弹性桩(或称柔性桩)之分，它们的受力特点是不同的，因而 应注意区分，即应先按下式计算出桩的变形系数：式中b0为桩的计算宽度 的单位为m-1当h2=1时，为弹 性桩。如图1所示，抗滑桩在滑坡推力和土体抗力作用下产生弹性挠曲变形，其受荷段.A) 的挠曲微分方程为采用等量分段h由桩顶往下至滑动面将桩长离散化，如图2(a)(b)所示。则任意节点i处的控制差分方程为斗 i- - , + &，, + ）：. + 提 I M 职烦=6 mlhEl普一 I 值 M + （7, -V.料It l由于桩顶自由，故在桩顶处弯矩和剪力均为0,即由

10、QJ得-r, =也与二o+2孔F +为二o及由 土x得J -2 二 2 I - 2% + 当5）中i=0,并将式（4）和式（5）代入可得（3），令式令则式（4）和（6）可以表达为r-1 = -iTo -6-iTin 二X + 4仿此，可得N =色j+乌N）Z-!=色一一 1 +弓 妃 Ji t q将（8）（9）代入（3）整理可得上式中令rh4-构一 i + 6+ 玲%-也 | +4 /d . b. = X/de =能一 J +L 一仇曲 I -则式（10 ）仍可写成式（7）的形式。滑动面以下，即锚固段可视为桩顶受水平荷载（集中 荷载剪力和弯矩）的侧荷桩，其挠曲微分方程为EJ乌 + 5 = 01

11、1）从桩底由下往上至滑动面也采用等量分段h见图2（c），则可得该段任意节点j处的控制差 分方程为?成_-如T十成一4小广。桩底按自由端考虑，由弯矩和剪力均为0的边界条件，同理可得到式（4）和式（5）。同样地式（12）中令j=0，可得将式（4）和（5）代入上式，整理可得令则防上2式可得将（14） （15）代入（12）整理可得（15），可得到下面的线性方程组由于抗滑桩在滑动面处应满足位移、转角、弯矩和剪力的连续条件。参照图2（b），（c），并 仿照式（7）（9）和式（13）解此线性方程组，可求得滑动面处及附近节点（上、下段各5个）处的位移。从而采用迭代 法可求得桩身各节点处的位移及其内力。计算结果

12、及准确度检验惜筮方海与有限莲分凭计算靖果比攻Table 1 Cdcuhtion ieuIu comparison between conventional method and finite difitrEtice method力法滑百处狩底姓 fcft/jll滑面妙弯雄/ Wlii涓面处W力八X最大有艳 / kNU.财力/ LX耕底处身甲 /k*ui拼底处 此W kN佐统中法110160.C03 3-0_0fl213 014.B16 66L5R13.202 105.911 5 .XI0.-9H020茬讣法OJJ15 90.0033 d-0JJ019H1906一知3 015.696 664.4513252 10(1 SL5.i000计算结果几乎完全一致，说明本文方法是正确可靠的。传统的查表手算法在桩底处出现了与 边界条件不完全吻合的残余误差，这是由于传统方法所用表中系数截断误差所致，而差分法 则没有出现这样的问题。6结论可以看出在实际问题中利用有限差分法对实际问题分析求解，再利用相应计算机程序如 FLAC、UDEC/3DEC和PFC程序，有时也可以根据自己需要根据算法编制