有限差分法及其应用

有限差分法及其应用

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,采用电壁(electric wall)边界条件,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处的控制差分方程为

即

本文来源：https://www.bwwdw.com/article/lyft.html