邓子平

导读

本文主要写给从事CAD/CAE/EDA/CFD等相关的软件研发测试人员，也可以作为有限元入门文章，对于从事仿真分析的工程师也可以对软件实现细节有所了解。

-----前言-----

从《越狱》里的胡克定理说起

看过美剧《越狱》的朋友对其中一个桥段应该印象很深刻，男主让狱友将事先画好的一张图铺在墙上，然后用打蛋器在固定点打孔，最后将墙敲穿推到。狱友问这是什么原理，男主说：这是利用了胡克定律。

高中物理课本对胡克定律的定义是：

弹簧在发生弹性形变时,弹簧的弹力F和弹簧的伸长量(或压缩量)x成正比。

很显然，如果仅仅用胡克定律无法解释几个小孔就能把墙凿穿。事实上男主角应该是事先拿到了墙的结构图，算好了砖块的距离，在砖头的缝隙处打孔。至于说原理的话，一张白纸很难从中间撕开 ，但如果纸中间有破损，沿着破损处就很容易撕开，非要扯上一个定律原理的话，孔边应力集中原理可以很好地解释这种现象。所以男主所说的胡克定律只不过是搪塞狱友的一个借口~

广义胡克定律

在材料的线弹性范围内，固体的单向拉伸变形与所受的外力成正比。

高中物理课本对胡克定律的定义可视为广义胡克定律的简化。

-----从刚体到非刚体-----

一根柔软的钢尺一端固定，一端悬空，悬空的一端会发生明显变形。但是高中物理只能告诉我们固定的一端受到向上的力以及力矩，无法告诉悬空一端的变形大小，这是因为其研究对象的属于刚体范畴，即研究的对象本身不会发生变形，而现实中绝大部分物体受力时自身会发生变形，工程力学中使用应力，应变，位移来度量对象的受力情况，而材料需要引入弹性模量和泊松比。有限元分析的几何对象是物体本身。

-----从一维到二维-----

弹簧是典型的一维单元，扩展到二维需要引入偏微分方程PDE。拉普拉斯方程和泊松方程是最简单的二维偏微分方程。这两个方程可以描述静电场，导热，二维波动，以及天体运动等。简单讲，PDE（Partial Differential Equation）是描述这个世界运行规律的一种方式。从量子原子微观运动，到宏观热，声音，电磁，以及天体万有引力规律，都有各自的方程来描述。声场对应赫姆霍兹方程，电磁波对应麦克斯韦方程，流体对应纳维斯托克方程，热传导对应泊松方程，量子力学对应薛定谔方程等等，而这个方程基本都是PDE。我们所学的各种自然定律规律，比如牛顿三定律，浮力定律等只不过是常微分方程的简化，而常微分方程则是PDE的简化。

有限元方法求解的对象是PDE，有限元方法的本质是将求解对象离散成小的单元，在每个小的单元上用基函数和形函数表征其坐标和物理量，使函数在其自由度所在的网格单元几何上满足偏微分方程的解（之所以是几何不是顶点，是因为自由度不一定在顶点）

-----有限“单元”篇网格-----

网格是个很广泛的称谓，为了避免歧义，FEM中的网格做如下规定：

1、网格的英文名为Mesh，非Grid

Grid 英文翻译也是网格，但是基本意思为格子，通常指结构化网格，即四四方方非常规整的网格。

2、网格是用于有限元计算的输入数据

网格对应于有限元的单元数据；线单元，面单元，实体单元都称之为网格（对，线单元的一条直线也是网格）

3、几何离散成的三角形称之为“面片”，非网格

大部分显示引擎底层都使用三角形来渲染对象，因此2维、3维几何都需要三角化（称之为“面片化”），比如一个长方体共有6个面，每个面需要离散成两个三角形，总共12个三角形。这12个三角形我们称之为12个“面片”，而非网格！一般情况下“面片”质量很差。“面片”有两个功能：一是用来做显示渲染数据；二是可以作为网格划分的输入数据，用来生成面网格。

4. 网格计算泛指“分布式计算”

“网格计算”（Grid Computing）是分布式计算的一种，它研究如何把一个需要非常巨大的计算能力才能解决的问题分成许多小的部分，然后把这些部分分配给许多计算机进行处理。这里的网格翻译也是Grid，和有限元计算无关！

-----单元-----

有限元方法的输入基本对象是“单元”（element)，即几何模型被离散之后形成的网格模型，通常根据单元的几何特点可以分为0维（点单元），1维（线单元）， 2维（面单元）和3维（实体单元）

实体单元（3维单元）很容易理解，因为任何物体都是3维的，非实体单元通常是对模型的抽象和简化。这种抽象和简化基于一定的前提假设，能在不降低计算准确度的前提下大幅降低计算时间。比如常见的梁单元，可以将实体梁简化为线单元，大大减少计算量。

弹簧只能在一个方向上发生变形，是典型的1维单元；同理壳单元（shell）需要XY两个方向来定义，是2维单元；四面体，六面体是3维单元，也称为实体单元；对象可以看做质点的为0维单元，比如称之为“定楼神球"的调谐质量阻尼器。

单元的阶数：

以三角形为例

一次线性单元为3个顶点

二次单元为6个点，包括3个顶点和3条边的中点

三次单元有9个点，每条边上有2个点

为了方便，统一将线性单元称之为0阶单元，二次单元称之为1阶单元，以此类推。

-----边界条件（Boundary Condition）-----

有限元求解的对象是偏微分方程（Partial Differential Equation），理论上偏微分方程的通解有无数多个，但实际上确定的工况下模型只有唯一一个确定解。而决定最后唯一解的就是边界条件！简单理解就是通解是一个带未知数的函数，而边界条件可以求解出这些未知数！未知数的个数和边界条件能确定的数值个数相同。

边界条件按照数据类型分为三类：

1、只有确定数值的，比如确定的电流，压力，温度，位移等数值，称之为第一类边界条件，英文为Dirichlet

2、用函数表示的，函数可以是用导数，积分表示的任意函数。比如随时间变化的荷载，电压；用导数表示的换热系数，称之为第二类边界条件，英文为Neumann

3、第一类和第二类的混合，称之为第三类边界条件，英文为Robin

很多书将三类边界条件用中文表示，但翻译有偏差，比如Dirichlet，翻译有狄里克雷，狄里克莱，狄立克莱，狄力克莱，狄里赫利。

对于工程中使用的英文名称，建议统一用英文或者无歧义的中文（第N类边界条件）表示，不使用英译。

边界条件的设置基于对有限元模型的正确理解。三维电磁分析时，设置的理想吸收边界需要将求解域包围住；施加的Port激励根据类型设置在不同的几何上；结构中的位移约束和荷载不能同时施加在同一几何上；流体中的进口速度压力和出口速度压力边界要与计算使用的CFD模型类型匹配。

在实际使用仿真软件时，用户无需过多关注第几类边界条件，因为在软件中边界条件统一使用离散的方式输入，即用户在设置边界条件时，通常是选中某个几何（比如选中一个面），然后通过对话框输入数据，将属性数值直接关联在该几何上，对于用函数表达的第二类，第三类边界，通常也是采用用户输入数据，采用简单拟合插值的方法得到表达式。同时商业软件中也会根据实际情况，提供多种便利的边界设置方法。比如要在100片相同的叶片上设置相同的荷载，可以先定义一个荷载，然后一起关联到所有几何上。这也是商业软件的价值所在：提高前处理效率！

在开发边界条件处理功能时，本质上是有共性的：即将数值和几何关联（最终是将数值和网格自由度关联），不同的只是数值的单位和类型。所以基类即可以实现大部分功能，节省开发工作！拟合插值功能也可以做成公用模块，插值方法无非线性，二次，对数插值常用的几种。

正确的使用单元和设置边界条件是得出正确计算的重要条件！

-----单元设计（针对软件开发）-----

ANSYS和ABAQUS中都有上百种单元，根据前面介绍，单元分类按照空间维度有0维，1维，2维，3维

按照面向对象设计方法，可以设计出四个接口基类；按照阶数来分，至少有0阶，1阶和2阶，阶数可以作为属性放在实现类中；引入自由度，单元设计会稍微复杂一点，每个节点自由度根据特定单元类型不同而不同，比如有些节点只考虑平动，不考虑转动，有些节点只考虑磁场不考虑电场；再引入多物理场，单元设计会更加复杂，以三维四面体单元为例，对于温度，0阶单元可作为默认单元，对于结构0阶单元误差过大，默认为1阶单元，而对于电磁分析由于自由度在边不在点上，0阶1阶2阶都适用，可根据几何实际情况选择。

早期的有限元分析程序多采用FORTRAN编写，没有采用面向对象思想和方法，面向过程的好处是程序流程会比较清晰，任何时候都能对单元的类型结构状态数据有准确的定位，但是难以维护和扩展。当使用类似面向对象C++语言设计大型有限元程序时候，不能抛弃面向过程的设计方法！过多依赖面向对象，从工程角度讲会造成软件开发的不可控，虽然从程序角度看会带来程序的灵活性，扩展性！

对于大型通用有限元分析程序，合理的设计单元架构是有效建立有限元模型的关键，单元框架需要同时供前处理和求解器使用！

-----误差-----

误差天然存在，误差不是错误，但是会导致错误

以实体建模为例， 网格划分精度通常为双精度，也就是 double型；而渲染显示（比如OpenGL）精度为单精度为float型。创建一个长方体，长方体的顶点坐标如果用鼠标从屏幕上取渲染数据值，即为float型数据而不是double型，生成的几何数据精度会丢失，直接结果就是生成的长方体会和其它相接触的实体产生干涉，如果不进行修复，划分网格就会在边界处生成非常细长的三角形，最终导致畸形网格，求解器无法进行计算。

有限元方法中有三种误差：

离散误差，累积误差和截断误差。

对于开发人员而言，需要注意累计误差和截断误差，比如等参单元使用的高斯积分，其积分误差和选择的积分点数和位置有关。在求解器求解线性方程组时，采用的不同的求解方法，其最终产生的误差也不相同。

网格的数目，密度，形状，不同的离散方法，不同的求解方法等最终产生的误差很可能导致完全错误的结果。对于研发而言，每一步产生的误差都需要仔细分析，给出处理逻辑。

-----求解器干了什么-----

有限元模型生成后，通常会导出保存为一个文件，文件中保存了节点，单元，材料，边界，荷载，工况，求解等设置信息，求解器会读入该文件进行计算。ANSYS为cdb，ABAQUS为inp，PATRAN为bdf。

一般商业软件都会把求解器做成单独的可执行程序（*.exe程序），单独启动进程求解，而不是集成在前处理器中。一是单独的进程便于管理和分布式计算；二是解耦前处理器和求解器，方便调试；三是可以方便导入到第三方工具中，缺点是模型文件容易被非法查看，可以通过二进制和加密来解决。

求解器读入有限元模型后，首先会检查模型的有效性，比如网格节点是否正确（重复节点，重复网格，索引错误等），还会检查网格质量，最小角度，边长，Aspect ratio，skew ratio 以及对应的求解配置是否正确，比如边界是否设置，工况是否合理（振动频率阶次过多），材料是否有效（比如泊松比大于0.5）等等。当出现不合理设置时，给出警告信息；出现错时，则返回不计算。模型检查合格之后，根据求解器设置，计算每个单元的刚度矩阵，然后组装刚度矩阵成全局矩阵，最终建立起一个非齐次线性方程组，最后求解该线程方程组。

在整个过程中最耗时的是求解方程组。受限于软硬件，求解大规模线性方程组仍然是世界性技术难题，目前对于上了十亿自由度的方程组，除了硬件GPU，分布式计算，超级计算机硬件加速外，并无快速有效的求解方法。显式计算方法不需要求解大规模线性方程组，更适合分布式计算。期待不久的将来量子计算能解决此类问题！

求解器计算完成后，会将计算自由度所在点的原始数值输出到文件。之所以说是原始数值，是因为求解器得到的可能并不是用户最终关心的数据，还需要进行一系列处理，比如有些求解器算出点的结果数据为积分点上的值，而非顶点值；用户通常关心Vonmises应力而非最大最小主应力；电路设计用户只关心SNP参数；天线设计用户更在意电磁场场强分布，等值面通常比单个点值更直观等等，所以需要对原始结果数据进行加工处理，给出对用户最有价值的结果！