微分方程数值计算方法简介

求解微分方程的常用数值方法有有限差分法、有限元法、有限有限体积法、边界元法、离散元法、光滑粒子流体动力学法等。

1.有限差分法

有限差分法(FDM)是一种直接将微分问题化为代数问题的数值方法。有限差分法具有简单、灵活以及通用性强的特点，容易在计算机上实现，是发展比较早且比较成熟的数值方法。

有限差分法的基本思想是先把问题的定义域进行网格划分，然后在网格点上，按适当的数值微分公式把定解问题中的微商换成差商，从而把表示变量连续变化关系的偏微分方程离散成为有限个线性代数方程，进而求出数值解。此外，还要研究差分格式的解的存在性和唯一性、解的求法、解的数值稳定性、差分格式的解与原定解问题的真解的误差估计、当网格尺寸趋于零时差分格式的解是否趋于真解等。

目前主要采用泰勒级数展开方法来构造差分格式,有一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等,其中前两种差分格式为一阶计算精度,后两种格式为二阶计算精度。

有限差分法已经成为求解各类数学物理问题的主要数值方法之一。在固体力学中,有限元法出现以前,主要采用差分方法; 在流体力学中,差分方法仍是较为常用的数值方法。代表性的商业软件有爆炸冲击分析软件 AUTODYN、岩土工程专业分析软件 FLAC。有限差分法适于求解边界较为规则的问题,如果边界条件比较复杂,就不太适用了。

2.有限元法

有限元法(FEM)是一种求解偏微分方程边值问题近似解的数值技术,求解时对整个问题域进行分解,每个子区域都成为简单的部分。

有限元法通过变分法,使得误差函数达到最小值并产生稳定解。其基本思想是将连续的求解区域分解成一组有限个小的互联单元,对每一单元假定一个合适的(较简单的)近解,然后推导求解这个域总的满足条件,从而得到问题的解。有限元法的实质是最小势能原理的应用,对每一有限元假定一个合适的、较简单的近似函数来分片地表示全求解域上待求的未知函数,单元内的近似函数通常由未知场函数或其导数在各个节点的数值及其插值函数来表达。如此一来,在一个问题的有限元分析中,未知场函数或其导数在各个节点上的数值就成为新的未知量(即自由度),从而使一个连续的无限自由度问题变成离散的有限自由度问题。一经求出这些未知量,就可以通过插值函数计算出各个单元内场函数的近似解,从而得到整个求解域上的近似解,达到用较简单的问题代替复杂的实际问题的目的。

有限元法的解不是准确解,而是近似解,因为实际问题被较简单的问题所代替。而且.随着单元数目的增加,即单元尺寸的减小,或随着单元自由度的增加及插值函数精度的提高,解的近似程度将不断改进。如果单元是满足收敛要求的,近似解将收敛于精确解。由于大多数实际问题难以得到准确解,而有限元法不仅计算精度高,而且能适应各种复杂形状因此有限元法成为一种广泛应用、行之有效的数值分析手段。

有限元法的缺点是难以应用于具有极大单元变形的情况,求解大型问题时需要的内存和计算量比其他数值算法要大得多、预估计算产生的误差比较困难。

采用有限元法编制的软件非常多,代表性的商业软件有 ANSYS Mechanical、LS-DYNA、MSC.Nastran 等。

3. 有限体积法

有限体积法(finite volume method,FVM)又称有限容积法、控制体积法，是计算流体力学中常用的一种数值计算方法,有限体积法基于积分形式的守恒方程而不是微分方程,该积分形式的守恒方程描述的是计算网格定义的每个控制体。有限体积法着重从物理观点来构造离散方程,每一个离散方程都是较小体积上某种物理量守恒的表达式,推导过程物理概念清晰,离散方程系数具有一定的物理意义,并可保证离散方程具有守恒特性。

有限体积法的基本思路是:将计算区域划分为一系列不重叠的控制体积,并使每个网格点周围有一个控制体积,将待解的微分方程对每个控制体积积分,便得出一组离散方程其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分,必须假定数值在网格点之间的变化规律，即假设值的分段的分布剖面。

有限体积法的优点是:

(1) 具有很好的守恒性

(2) 更加灵活的假设,可以克服泰勒展开离散的缺点。

(3)对网格的适应性好,可以很好地解决复杂的工程问题

(4)在进行流固耦合分析时,能够和有限元法完美融合。

大多数计算流体力学软件如 Fluent,STAR-CD、Flotherm 等采用的都是有限体积法此外有限体积法也是爆炸冲击计算软件 MSC.Dytran 的基本算法之一。

4 .边界元法

边界元法(boundary element method,BEM)是继有限差分法有限元法之后发展起来的又一种重要的数值方法。边界元法只在定义域的边界上划分单元，用满足控制方程的函数去逼近边界条件。与有限元法相比,边界元法降低了问题的维数和自由度数,具有单元数量少、数据准备简单、计算时间少等优点,已经成为现代科学和工程数值分析的有效工具。边界元法的基本思想是应用格林定理等,将问题的控制方程转换成边界上的积分方程然后引人位于边界上的有限个单元将积分方程离散求解。借鉴有限元法划分单元的离散技术,通过对表面边界进行离散,得到边界单元，将边界积分方程离散成线性方程、经过离散后的方程组只含有边界上的节点未知量,因而降低了问题的维数,最后求解方程的阶数降低,数据准备方便，计算时间缩短。另外,边界元法通过引用问题的基本解而具有解析与离散相结合的特点，使得计算精度较高。由于积分方程可以用加权余量法得到,这就避免了寻找泛函的麻烦。

边界元法的主要缺点是它的应用范围以存在相应微分算子的基本解为前提,对于非均匀介质等问题难以应用,故其适用范围远不如有限元法广泛,而且通常由它建立的求解代数方程组的系数矩阵具有非对称、稠密,甚至在某些情况下病态的特性。应用边界元法分析大规模问题,系统方程求解过程消耗大量计算资源和时间,对解题规模产生较大限制。对一般的非线性问题,由于在方程中会出现域内积分项,从而部分抵消了边界元法只离散边界的优点。

边界元法代表性的商业软件有三维边界元流体力学分析软件 ANSYS LINFLOW等在 LS-DYNA 软件中也包含了边界元声学计算功能。

5. 离散元法

离散法(discrete element method 或者 distinct element method,DEM)是 20 世纪70 年代由 Cundall 首先提出来的,起源于分子动力学,是为研究岩体等非连续介质的力学行为而发展起来的一种数值方法。

离散元法的基本原理是牛顿第二定律,其基本思想是把求解域分离为离散单元的set使每个离散单元满足牛顿第二定律,离散单元本身一般为刚体,单元间的相对位移等变形行为一般由联接于节点间的变形元件来实现。离散单元可以平移、转动或变形。各个离散单元在外界的干扰下就会产生力和力矩的作用,由牛顿第二定律可以得到各个离散单元的加速度，然后对时间进行积分,就可以依次求出离散单元的速度、位移,最后得到离散单元的变形量。离散单元在位移向量的方向会发生调整,这样又会产生力和力矩的作用。如此循环直到所有离散单元达到一种平衡状态或者处于某种运动状态之下,使各个离散单元满足运动方程,用时间步长迭代的方法求解各个离散单元的运动方程,继而求得不连续体的整体运动形态。离散单元法允许单元间的相对运动,不一定要满足位移连续和变形协调条件，计算速度快,所需存储空间小,尤其适合求解大位移和非线性的问题。

离散元法比较适合于模拟节理系统或者离散颗粒组合体在准静态或动态下的变形过程,如颗粒材料的混合、分离、注装、仓储和运输过程,以及连续体结构在准静态或动态条件下的变形及破坏过程。

目前,离散元的理论体系并不完善,其局限性表现在接触模型和参数难以确定、稳定性难以保证等方面。

基于离散元法的商业软件有 UDEC、PFC、GRANULE、EDEM 等,离散元也是 LSDYNA 软件的重要功能之一。

6 .光滑粒子流体动力学法

光滑粒子流体动力学(smoothed particle hydrodynamics,SPH)法是一种为解决天体物理中涉及的流体质团在三维空间无边界情况下的任意流动问题的纯拉格朗日方法。后来发现解决如连续体结构的解体、碎裂、固体的层裂、脆性断裂等物理问题,以及产生大变形的流体动力学问题、考虑材料强度的固体动力学问题也非常有效。总之,SPH 法适合用来求解具有大变形、复杂边界和物质交界面等复杂的单相或多相流体动力学问题,是最早的无网格方法。

光滑粒子流体动力学法的理论基础是插值理论,通过一种称之为核函数的积分核进行插值近似.从而将流体力学方程转换为数值计算用的 SPH方程组。光滑粒子流体动力学法的基本思想是用一系列的粒子来表示求解域,这些粒子具有流体的密度、速度、热能等物理量,粒子之间不需要任何连接,即具有无网格特性,粒子的密度、位移、速度、压力等物理量的更新只与时间有关:用积分近似和粒子近似离散流体动力学控制方程,产生带状或离散化的稀疏系数矩阵,某一粒子场函数的值通过支持域内相邻粒子的叠加求和计算得到;由于所使用的用于求和的局部粒子为当前时间步的粒子,所以 SPH 法具有较强的适应性;将粒子近似法应用于所有偏徽分方程的场函数相关项中,可得到一系列只与时间相关的离散化形式的富徽分方程,应用显式积分法来求解常微分方程以获得最快的时间积分,并可得到所有粒子的场变量随时间的变化值。

SPH法相对于传统的基于网格的方法具有一些特别的优点:具有自适应性、无网格性以及拉格朗日公式与粒子近似法的结合。一方面,由于 SPH 法的近似过程不受网格的限制，因此可避免网格极度大变形造成的精度下降问题;另一方面，SPH 法是一种拉格朗日方法,每个粒子点代表一种独立的物质,可以自然且直观地跟踪粒子所代表的物质属性,能够方便地指述多物质流动问题。

SPH法也存在一些数值方面的困难,如拉伸不稳定性会引起数值断裂、形函数不一致性、引入本质边界条件存在困难,需要比较复杂的接触算法等。

基于SPH法的商业分析软件有美国 CentroidLab 公司开发的 NEUTRINO,开源 SPH软件有 DualSPHysics,SPHinXsys、SPHERAL 等,在 LS-DYNA、AUTODYN、ABAQUS软件中均包含 SPH计算功能。

摘自：《Ls-dyna有限元建模、分析和优化设计》辛春亮