有限差分法及其应用

有限差分法及其应用1有限差分法简介有限差分法（FDM）是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方程将解域划分为差分网格，用有限个网络节点代替连续的求解域。有限差分法通过泰勒级数展开等方法，把控制方程中的导数用网格节点上的函数值得差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。2有限差分法的数学基础有限差分法的数学基础是用差分代替微分，用差商代替微商而用差商代替微商的意义是用函数在某区域内的平均变化率来代替函数的真是变化率。而根据泰勒级数展开可以看出，用差商代替微商必然会带来阶段误差，相应的用差分方程代替微分方程也会带来误差，因此，在应用有限差分法进行计算的时候，必须注意差分方程的形式，建立方法及由此产生的误差。3有限差分解题基本步骤有限差分法的主要解题步骤如下：1）建立微分方程根据问题的性质选择计算区域，建立微分方程式，写出初始条件和边界条件。2）构建差分格式首先对求解域进行离散化，确定计算节点，选择网格布局，差分形式和步长；然后以有限差分代替无线微分，以差商代替微商，以差分方程代替微分方程及边界条件。3）求解差分方程差分方程通常是一组数量较多的线性代数方程，其求解方法主要包括两种:精确法和近似法。其中精确法又称直接发，主要包括矩阵法，高斯消元法及主元素消元法等；近似法又称间接法，以迭代法为主，主要包括直接迭代法，间接迭代法以及超松弛迭代法。4）精度分析和检验对所得到的数值进行精度与收敛性分析和检验。4商用有限差分软件简介商用有限差分软件主要包括FLAC、UDEC/3DEC和PFC程序，其中，FLAC是一个基于显式有限差分法的连续介质程序，主要用来进行土质、岩石和其他材料的三维结构受力特性模拟和塑性流动分析；UDEC/3DEC是针对岩体不连续问题开发，用于模拟非连续介质在静，动态载荷作用下的反应；PFC是利用显式差分算法和离散元理论开发的微、细观力学程序，它是从介质的基本粒子结构的角度考虑介质的基本力学特性，并认为给定介质在不同应力条件下的基本特征主要取决于粒子之间接粗状态的变化，适用于研究粒状集合体的破裂和破裂发展问题，以及颗粒的流动(大位移）问题。5有限差分法的应用1）光子晶体光纤色散的有限差分法研究光子晶体光纤的横截面结构和特征参数如图1所示。对结构如此复杂的光子晶体光纤,采用有限差分法加以研究。有限差分法是光波导分析中广为采用的数值方法之一,它利用台劳展开用有限差分式代替波动方程中的微分式,得到关于场分量的有限差分方程。在纵向(z方向)均匀分布的光子晶体光纤中,电场的横向分量Ex和Ey满足如下耦合方程组其中,Ex和Ey为电场的两个偏振分量,εr为相对介电常数,k0为真空中的波数,β=koneff为传播常数,neff为模式有效折射率。在光子晶体光纤分析中,通常求解基模(类比于传统光纤,也称为HE11模)的两个偏振分量。文献中已经证明了横向结构对称分布的光子晶体光纤基模的线偏振特性,因而(1),(2)式中的电场偏振分量Ex和Ey可以认为不再耦合,此时(1),(2)式将变为两个独立的分别关于两个偏振分量Ex和Ey的半矢量形式的波动方程(3),(4)式中保留左边的第二项对于研究空气填充率较高的光子晶体光纤是必要的,因为这时标量波动方程给出的解会产生较大的误差。在x,y方向采用相同的网格间距h,利用五点差分格式可以建立基于半矢量波动方程(3),(4)的差分方程。其中E为Ex或Ey,p,q为网格点编号。对于Ex,各项系数分别为对于Ey,各项系数分别为(5)式可化为特征值方程其中A为系数矩阵,E为由各网格点构成的向量。计算结果为将基于半矢量波动方程的差分法得到的结果与文献中全矢量方法得到的结果比较,取与文献中相同的求解参数:光纤参数为空气孔间距Λ=2.3μm,空气孔直径d=1.0μm,整个求解区域取为6Λ×6Λ,石英的折射率为1.45,波长为1.5μm,采用电壁(electricwall)边界条件,x,y方向格点数均为240。基于半矢量波动方程的差分法得到结果为neff=1.42806,与文献中给出的结果1.42868以及其他方法得到的结果较为接近。考虑到全矢量方法计算结果表明在求解某一偏振态时另一偏振分量并不为零,但相差几个数量级,而半矢量波动方程完全忽略了另一偏振分量的影响,以及网格剖分的方法不同,这一结果是较为精确的。2）抗滑桩全桩内力计算“m-k”法的有限差分法如图1所示，设抗滑桩全长为H，其中滑动而以上，即受荷段桩长为h1;滑动而以下，即锚固段桩长为h2。设滑坡推力ET按梯形分布，地而和滑动而处的分布荷载集度分别为q0和qa。桩前滑体的剩余抗滑力为E'。桩的截而尺寸:长x宽为axb，桩的抗弯刚度为EI。为便于公式的推导，设桩前滑体地基系数为K1(z)=mZ，滑而以下锚固体地基系数为K2(z)=K(常数)。由于抗滑桩有刚性桩和弹性桩（或称柔性桩）之分，它们的受力特点是不同的，因而应注意区分，即应先按下式计算出桩的变形系数：式中b0为桩的计算宽度的单位为m^-1当h2=1时，为刚性桩;h2=1时，为弹性桩。如图1所示，抗滑桩在滑坡推力和土体抗力作用下产生弹性挠曲变形，其受荷段.A)的挠曲微分方程为采用等量分段h由桩顶往下至滑动面将桩长离散化，如图2(a)(b)所示。则任意节点i处的控制差分方程为即由于桩顶自由，故在桩顶处弯矩和剪力均为0，即由得及由得令式（3）中i=0，并将式（4）和式（5）代入，可得令则式(4)和（6）可以表达为仿此，可得将（8）(9）代入（3）整理可得上式中令则式（10）仍可写成式（7）的形式。滑动面以下，即锚固段可视为桩顶受水平荷载（集中荷载———剪力和弯矩）的侧荷桩，其挠曲微分方程为从桩底由下往上至滑动面也采用等量分段h见图2（c），则可得该段任意节点j处的控制差分方程为即桩底按自由端考虑，由弯矩和剪力均为0的边界条件，同理可得到式（4）和式（5）。同样地式(12)中令j=0，可得将式（4）和(5)代入上式，整理可得令则防上2式可得将（14）（15）代入（12）整理可得由于抗滑桩在滑动面处应满足位移、转角、弯矩和剪力的连续条件。参照图2（b），（c），并仿照式（7）~（9）和式（13）~（15），可得到下面的线性方程组解此线性方程组，可求得滑动面处及附近节点（上、下段各5个）处的位移。从而采用迭代法可求得桩身各节点处的位移及其内力。计算结果及准确度检验计算结果几乎完全一致，说明本文方法是正确可靠的。传统的查表手算法在桩底处出现了与边界条件不完全吻合的残余误差，这是由于传统方法所用表中系数截断误差所致，而差分法则没有出现这样的问题。6结论可以看出在实际问题中利用有限差分法对实际问题分析求解，再利用相应计算机程序如FLAC、UDEC/3DEC和PFC程序，有时也可以根据自己需要根据算法编制一个计算程序求出方程组的解从而得出的结果与查表结果和实际结果对比，精确对很高，从而证明了其正确性与实用性。7参考文献：1《计算机在材料科学与工程中的应用》-中南大学出版社、2《抗滑桩全桩内力计算“m-k”法的有限差分法》戴自航，彭振斌3《光子晶体光纤色散的有限差分法研究》栗岩锋,刘博文,王子涵,胡明列,王专,王清月(天津大学精密仪器与光电子工程学院光电信息技术科学教育部重点实验室,天津300072)