结构力学分析方法的江湖

weixin

现阶段结构力学分析的江湖几乎已被有限单元法一统江山。但实际上，有限单元法只是数值方法求解微分方程的众多手法之一，在此文中，将介绍现阶段几位一样能够参与结构分析华山论剑的绝顶高手。

武林盟主 —— 有限单元法

如今江湖，各路豪杰并起，热闹非凡，上至通用分析软件中鼎鼎大名的ABAQUS、ANSYS、SAP2000等，下至结构工程师的老朋友PKPM、YJK、Midas，若要追根溯源，有限单元法是他们共同的基础。在结构设计界，据我所知，有限单元法已经一统江湖，号令天下，莫敢不从，俨然武林盟主之态。那么，这位盟主出身如何呢，且看以下分解。

盟主的身世

有限单元法的发展

有限单元法思想最早出现在1943年，Courant 采用定义在一系列三角形区域的分片连续函数结合最小势能原理求解St.Venant 扭转问题。由于当时计算工具的限制，这种方法并没有在接下来的时间里得到较好的发展。直到1960年以后，电子计算机的广泛应用和发展，现在意义上的有限单元法才真正发展起来。

现代意义上的有限单元法第一个成功的尝试出现在1956年，这是由Tuner 和Clough 等人在分析飞机结构时提出的。他们使用有限单元法三角形单元，得出了平面应力问题的正确解答。1960年Clough 将其进一步推广到平面弹性问题，并第一次提出“有限单元法”这一名称。而在这之后的1965年，O. C. Zienkiewicz 和Y. K. Cheung 发现存在变分形式的所有场问题即满足变分原理的所有场问题，都可以用于固体力学有限元法相同的步骤求解。

但是并不是所有问题都能找到对应的变分原理，在不能找到变分原理的问题中应用有限单元法，就必须找到能够给出有限元格式的方法。在1969年，Oden J T 从所谓的能量平衡法出发，得出热弹性问题的有限元分析方程组。1969年，B. A. Szabo 和G. C. Lee 将加权余量法，特别是Galerkin 法引入有限单元法导出标准的有限元过程来求解非结构问题。加权余量法引入有限单元法后，使得有限单元法的应用范围得到了极大拓展。而且有限单元法同门中还有两位骨骼惊奇的奇侠，一般的有限单元法，形函数都是一阶线性的，对于复杂的几何边界条件比较吃力，所以有了等参元法和Wilson 非协调元，这两位形函数中含二次项的选手。

20世纪70年代，有限单元法的发展已达到高峰，其基本理论和方法已趋成熟。

盟主的招式

有限单元法的基本步骤

盟主有限单元法在解决敌手时总有个固定的套路。

起手式，将求解域离散为一系列单元。就是现在各路软件中划分网格的过程。

紧接一招隔空猜物，就是在各个单元内用基于单元各结点值的单元插值函数来拟合此单元上的未知场函数，我们常见的力学求解问题，一般都是先求解整体结构的位移函数。

然后一招万法归宗，根据力学原理求出各单元的刚度矩阵。再由单元刚度矩阵凝聚为总刚矩阵，配合边界条件解出数值解。

在一般的结构力学问题中，这几招下来基本无坚不摧。

盟主的罩门

有限单元法的瑕疵

盟主虽然已然一统江湖，但是其余高手仍未被除尽，可见盟主武功中也存在弱点。

首先就是单元的划分，虽然许多有限元程序有了自动划分单元的功能，但至今还无法保证划分的网格一定合理，由于网格划分不合理导致计算精度差甚至不收敛也不能完全避免。

其次，现有的有限单元法大都采用一阶线性形函数，在拟合复杂问题时效率不高，而等参元法拟合的二次形函数并不完备，Wilson 非协调元由于不符合变形协调条件，不能保证收敛。

再者，有限单元法在涉及大变形问题时，由于计算中网格会出现严重的扭曲，这样不但需要网格的重构，而且也会导致计算精度的下降。在处理裂纹发展问题时，由于裂纹发展的不确定性，在计算中随着裂纹的发展不得不不断地重新划分单元来模拟整个裂纹发展的动态过程。

隐侠 —— 无网格法

隐侠的由来

无网格法的发展

针对盟主的罩门，另一位高手修炼了另一门绝学——无网格法，其计算与网格无关。这种无网格法的出现，在大变形和裂缝发展问题中，不像有限单元法一样需要重新划分网格，因此，无网格法在处理这类问题时不仅能够提高计算精度，并且还能减小计算难度。

无网格法的出现最早可以追溯到20世纪70年代。当时人们开始研究对非规则网格有限元差分法的研究，这是无网格法的前身。1977年，Lucy 提出了“Smoothed Particle Hydrodynamics”方法，此方法即被简称为SPH方法。1995年，Liu 等人在其基础上提出了再生核粒子方法——“Reproducing Kernel Particle Method”简称为RKPM 方法。RKPM 方法和SPH 方法一样都是基于核近似的无网格方法，它在SPH 方法基础上引入修正函数施加再生条件并且采用高斯积分。RKPM 方法不仅解决了SPH 方法在边界条件上不一致性的问题，而且完全消除了SPH 方法的张力不稳定。但由于RKPM 方法不满足Kronecker delta 函数性质，这使得这种方法施加边界条件比较困难。

除了基于核近似的无网格法以外，还有采用移动最小二乘法近似的移动最小二乘法——“Moving Least Square”简称为MLS方法。在1992年，Nayroles 等人将Galerkin 法引入移动最小二乘法，从而提出了弥散单元法——“Diffuse Element Method”简称DEM 方法。此后的1994年，Belytschko 等人对弥散单元法进行了改进，在形函数的导数中保留了被弥散单元法省去的项，并且引入Lagrange 乘子法应用边界条件，这就是无网格Galerkin 法——“The element-free Galerkin Method”简称为EFG 方法。

其后又有诸多如有限点法——“The Finite Point Method”简称FPM 方法、无网格局部Petrov-Galerkin 法——Meshless Local Petrov- Garlerkin Method，简称MLPG 方法、再生核阶谱单位分解方法等方法。

隐侠的招式

无网格法的基本思想

鄙人对无网格法涉猎不深，从基本思想上来讲无网格法是在求解域上布置有限数量的结点，在每个结点上定义覆盖结点周围一定区域的结点权函数或称结点核函数来拟合该结点附近某一区域内的物理量。每个结点所影响的区域称为其支撑域（支撑域可部分互相重合），由所有结点的支撑域来覆盖整个求解域。这样通过各个结点的权函数的叠加就可以拟合待求解函数。

为何归隐

无网格法的缺点

隐侠无网格法一身奇绝武艺，甚至摆脱了有限单元法网格的桎梏，可以无视网格划分是否合理达到无招胜有招境界，可是为何现实中却无一款商业软件采用此法呢？归根结底是由于隐侠招法过于奇诡，难以驾驭。就以我们熟悉的结构分析为例，我们已经习惯于在模型的一些的节点上愉快的施加边界条件。给一个节点刚接约束节点就不能有位移的畅快体验，在无网格法中几乎成为泡影。因为无网格法中节点的支撑域重合，通常不满足Kronecker 函数特性，即一个节点上的位移并不仅由此节点的形函数确定，周围的节点都可能会对其有影响，所以就算要施加一个位移为0的约束也是千难万难，而且无网格法的积分计算较有限元法复杂，所以这位高手暂时只能在特定领域内栖身，江湖中往往仅仅只留下他的传说。

盟主的觊觎者 —— 广义有限单元法

武林中从来不缺惊才绝艳的新人。广义有限单元法就是这么一位可能承继盟主之位的新人。广义有限单元法最早是在20世纪末由Babuska I 等人提出，将传统有限单元法中的插值函数作为单位分解函数，从中我们能够看出广义有限单元法实际上是传统有限单元法和无网格法中单位分解法结合产生的。广义有限单元法对比于传统有限单元法最大的区别就是，将能够反映待求解问题特性的函数引入到插值函数中，从而能避免网格重构，提高计算精度与计算效率。

在大致相同时期，我国也有学者提出了类似的思想。梁国平教授在1995年借鉴石根华的流形思想，提出了广义有限元方法。此后，栾茂田等人引入广义结点的概念，也就是在传统有限元方法的结点中增加新的自由度，从而得到不增加结点就可以模拟高阶变形的广义结点有限元方法。辛娅云在程尧舜教授的指导下改进了栾茂田等人提出的广义结点有限元方法，具体做法是将各个结点的局部坐标原点设置为该结点本身，给出了新的插值函数。在2007年，S. Rajendran 和B.R. Zhang 综合了最小二乘法、有限单元法和点插值法，创造了一种新的四结点四边形单元，并应用该单元解决了平面应力问题。之后，又将这种单元应用于自由振动分析和非线性柯西材料的平面应力分析。此后，李刚在程尧舜教授的指导下又在辛娅云的广义结点有限元方法上进一步改良，引入邻近点概念结合点插值法提出了半弥散单元法。半弥散单元法是在原有的有限单元法形函数上加入新的自由度使之成为二次完备的形函数，在拟合复杂问题时较传统有限单元法收敛的更快，精度更高。

初试身手

半弥散单元法的应用

本部分是半弥散单元法展示实力的部分，本节中的符号约定：

   · SDEM——表示本文中半弥散单元法；

   · T3——表示传统三结点三角形单元的有限单元法。

本节以数值算例比较半弥散单元法与传统三结点三角形单元，在稳态热传导问题计算过程中的优劣。

四边边界温度已知的矩形温度场问题
如下图所示一无限长矩形柱体的横断面，边长分别为L1 与L2，材料性质均匀，且柱体内部无热源。由于柱体无限长，而边界上热交换条件与纵向坐标无关。因此，问题简化为二维稳态热传导问题。

无限长柱体截面及边界条件

该问题的数学模型可以用以下微分方程表示：

边界条件为

区域内的温度分布解析解为：

为了方便计算，本文中将L1、L2 都取为1，边界条件中TW 也取作1；同时，材料的热传导系数Kx、Ky 均取为1；选取计算区域中心位置P(0.5,0.5) 为考察点。

根据问题的解析解式（上式）可以求出P 位置处温度值应为

半弥散单元法和传统有限单元法在计算中采用同样的网格划分。网格划分如下图所示。

网格Ⅰ(2*2)

网格Ⅱ(4*4)

网格Ⅲ(8*8)

网格Ⅳ(12*12)

采用以上四种网格划分方法，分别采用半弥散单元法和传统有限单元法求解该问题，求得的考察点P 温度值与用区域内的温度分布解析解求得的考察点P 温度值解析解列于下表。

表1 考察点P在两种方法不同网格划分下计算结果汇总

四种网格划分方式下相对误差对比可以下图表示。

四种网格划分下相对误差

如以自由度数作为比较标准。在表2中列出了半弥散单元法和传统有限单元法自由度数对比。

表2 两种算法各网格划分下的自由度数

从表2可知，半弥散单元法的自由度数与传统有限单元法的自由度数相同，这是由于将原广义结点有限元中增加的自由度可用结点信息表示，在形成刚度矩阵时凝聚掉。

下面进行基于自由度数目的误差分析，由于自由度数目在网格划分加密的过程中并不是线性增加的，在分析中我们将对自由度数目取自然对数作为衡量指标。

两种方法自由度数目与相对误差

可以看出半弥散单元法无论是在相同网格还是在相同自由度数目下，在计算的收敛速度和精度方面都远远优于传统有限单元法。

乱点英雄谱

观今时结构分析手段之江湖，豪杰并起，百舸争流。而有限单元法凭借其成熟的理论，严整的招法，为各大软件开发商所喜，因之一统江湖，一时之间风光无两。但其严整的招法既有好处也带来了缺陷，对于一些大变形问题，由于对手招式过于诡奇，应付起来难免措手不及。而无网格法招式轻灵多变，应付有限单元法所不能解决的对手时长袖善舞，应付自如，但也因自身套路难以琢磨而难以为众人所用。广义有限单元法则有些集二者之长的意味，但作为新晋高手，至今仍未显大名于世。

本人学识粗鄙，天下英雄实不能尽识，如有错漏，还望各位看官海涵。

参考文献：

[1] 王勖成.有限单元法.北京：清华大学出版社，2003

[2] 陈晓珞. 二次广义有限元性质及其边界条件的处理方法：[硕士学位论文]. 上海：同济大学，2011

[3] 李刚. 半弥散单元法：[硕士学位论文]. 上海：同济大学，2007

来源：非解构微信公众号（ID：non-structure），作者：无锋。