一篇文章入门多物理场有限元（全篇）

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

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

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

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

弹簧在发生弹性形变时,弹簧的弹力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参数；天线设计用户更在意电磁场场强分布，等值面通常比单个点值更直观等等，所以需要对原始结果数据进行加工处理，给出对用户最有价值的结果！

-----什么是多物理场-----

工程领域的多物理场通常是指宏观领域的物理场，包括电磁，流体，结构，热，声，光（光本质上属电磁）等等，偏向于工程应用领域。需要和理论物理中的四大物理场（微观领域的强力，弱力，以及电磁力和天体的万有引力）相区别。

-----形函数，基函数，伽辽金-----

伽辽金，俄罗斯人，伽辽金方法一如俄罗斯的武器，便宜，耐用，适应能力强，保养简单。

形函数（shape function）是用一个试探性的函数来表达单元内的偏微分方程，通常取多项式。

基函数（basis function）是用来表达一个单元内的物理量的坐标表达式。

伽辽金方法是加权余量法的一种（加权余量包括 配点法 子域法 矩量法 最小二乘法 伽辽金法），将形函数和基函数设为相同，同时将偏微分方程强形式变成弱形式，降低了求解限制，对于非连续介质，只要在插值点满足要求即可，提高了其适用范围。相对于其他方法，伽辽金方法最终生成的线性方程组的系数矩阵，稀疏对称，更容易求解。简单讲，伽辽金方法简化了有限元方法使用门槛，搭建起了有限元数值理论到软件实现的桥梁，成为早期商业软件实现的基础。

-----等参单元/高斯积分-----

等参单元：划分网格后，一般网格并不是标准网格（等边三角形，等边四面体或者立方体），直接使用形函数会造成计算上的困难，所以通常使用一个新标准网格来代替原有网格，并在标准网格和原有网格之间建立映射关系，新的网格称为等参单元。

高斯积分：在使用等参单元后，会有很多定积分需要求解，通常方法是在积分区域内按照一定规则选出一些点，称为积分点，算出被积函数在积分点的值，然后再乘以加权系数作为近似积分值。高斯积分是精度最高的一种方法。这也是为什么很多有限元软件（比如ABAQUS），最终结果并不在顶点上，而是在积分点上，需要注意转换。

--显示算法/隐式算法/稳态/瞬态/动力--

显示算法/隐式算法在比较结构动力分析方法中比较常见，显示算法中的节点状态主要由前两部决定，不需要求解整体刚度矩阵，适合大规模并行计算，ABAQUS Explicit和LSDYNA主要使用显式方法。

通常说的稳态指没有时间项的计算，物理场处于平衡不变化的状态。

瞬态/动力 中有时间项，即物理量随时间变化而变化，温度热传导，冲击荷载，结构振动，电磁时域都是和时间相关。常用求解方法是将时间离散，求解每一个时间点的状态。时间的离散步长需要根据实际情况计算，时间步长过大，求解不精确，过小，浪费计算资源。

几个多物理场例子

--被风吹塌的桥（流体，结构，振动）--

塔科马海峡吊桥（英语：Tacoma Narrows Bridge）是位于美国华盛顿州塔科马的两条悬索桥，横跨塔科马海峡，第一座桥于1940年建成，但不到五个月便倒塌，公认的直接原因是桥梁在受到强风的吹袭下引起卡门涡街，使桥身摆动，当卡门涡街的振动频率和吊桥自身的固有频率相同时，引起吊桥剧烈共振而崩塌，这次事件成为研究空气动力学卡门涡街引起建筑物共振破坏力的活教材，也被记载为20世纪最严重的工程设计错误之一。而事故的根本原因是节省成本，减少桥面厚度所引起的刚度降低，虽然桥梁整体强度毫无问题，但是刚度降低破坏了整体稳定，无形中增加了风载造成的冲击，为事故埋下了隐患；事实上即使桥梁不被吹垮，长期经受高强度的冲击荷载，其使用寿命也会大大降低！

卡门涡街是流体力学中重要的现象，在自然界中常可遇到，在一定条件下的定常流体绕过某些物体时，物体两侧会周期性地脱落出旋转方向相反、排列规则的双列线涡，经过非线性作用后，形成卡门涡街。如水流过桥墩，风吹过高塔、烟囱、电线等都会形成卡门涡街。常用的CFD软件中都会用卡门涡街作为基准测试建立标准算例。

--NVH（声音，振动，声振舒适度）--

NVH是衡量汽车舒适度的一个重要指标，在中国，随着家庭汽车越来越多，NVH（Noise、Vibration、Harshness）在汽车研发中成为热门研究内容，汽车行业的火爆可是说是NVH发展的重要推手。有统计资料显示，整车约有1/3的故障问题是和车辆的NVH问题有关系，而各大公司有近20%的研发费用消耗在解决车辆的NVH问题上。

-----任正非采访（热分析）-----

华为老总任正非在一次接受采访时表示散热和发热机理可能是电子技术最核心的竞争力

“芯片越来越小，元器件越来越集中，硬件工程、电子工艺最大的问题就是散热。有专家表示，未来50%的能源将消耗在芯片上，散热和发热机理也可能是电子技术最核心的竞争力”，未来电子产品不仅要进行电磁兼容，电源完整，信号完整，机械强度等仿真，散热分析更会成为重中之重。

-----1亿是1千万的10倍？-----

早期笔者使用矩量法求解线性方程组，在不使用快速多级方法，自由度达到3万时，台式机上已经无法求解出，8G的机器内存不够用。对于满秩矩阵的线性方程组，常规求解方法时间复杂度为n^3(n的3次方)。

对于自由度1千万以下稀疏矩阵的求解，好的台式机基本能应付，而当自由度达到1亿的时候，简单的将硬件乘以10倍完全不能满足要求。因为计算的空间复杂度，时间复杂度并不是线性，通常是NlogN，N^2或者更高。当自由度达到1亿时，不仅需要对硬件核心部件CPU，内存扩容，而且在磁盘阵列，I/O，并发计算，GPU，网络，带宽等方面都提出了更苛刻的要求。

在算法方面看，减少网格密度，在物理量梯度大的地方加密，无变化的地方将网格变稀疏，可以有效减少计算量；另外优化求解算法本身，使其更加易于并行化计算。

从目前来看，计算机硬件计算能力的更新速度，跟不上指数级求解规模的增加速度。这也是量子力学发展的最大推动力！

-----设计和优化-----

关于设计和优化，之前谈过很多。在工程项目中，有限元方法仿真只是其中一环，即对模型进行计算得到仿真结果。按照ANSYS的说法，多物理场仿真仅仅停留在第二个层次（第一层次是单物理场仿真）。ANSYS现在关注更高第三层次的设计驱动仿真 和 第四层次的业务驱动仿真。简单讲设计驱动仿真让仿真来决定设计，而业务驱动仿真则是让仿真决定业务，之前做过详细介绍，可以查看历史文章。

其实按照目前的水平，在世界范围内，能做到第三层次即仿真驱动设计，或者优化设计，已经是把仿真做的相当不错了！

--不同物理场的有限元方法一样吗--

理论基础一样，但实现方法不一样

有限元的核心之一在于对形函数的选择，也就是形函数要能准确表达PDE以及边界条件。以电磁单元为例，早期使用和结构相同的三角面片单元（自由度在定点），但是在计算边界的时候发现电磁场不满足连续条件，后来发明了矢量单元（自由度在边上）来解决这个问题，形函数和结构分析使用的形函数不同。

再比如热分析需要对环境温度进行定义，电磁分析需要设定吸收边界（物体到吸收边界区域也需要划分网格），而结构中则不需要。结构分析中常见的沙漏，自锁，电磁分析中因为自适应网加密数目过少造成的伪解，CFD边界层网格没有加密造成的求解误差过大，都是由各自不同物理场的特点决定的。

--多物理场软件--

COMSOL是名至实归的多物理场软件，在前处理方面提供了多种物理场建模功能，同时可以自定义偏微分方程，求解偏微分方程，这是其它软件无法比拟的的。COMSOL作为学习，科研绝对是利器！所以在高校推广的特别好，在研究所研究院也有相当的市场！But，在工程领域，COMSOL拓展的比较艰难。这也是本人经常说到的一个问题，软件只是工具，一旦涉及到工程领域，软件必须要能解决实际问题。从工程角度看，COMSOL就像是一个大杂烩，似乎什么都能做，但什么都做不精通！涉及电磁电子有专业的EDA工具，做结构有老牌的ANSYS，ABAQUS，空气动力飞行器有老牌的NASTRAN，流体Fluent，Star-CD，材料化学有VASP，MS。

综合来看COMSOL和MATLAB属于同一类型的产品，着重在于研究解决某个点上的问题，而非解决综合工程问题

----ANSYS----

作为老牌商业仿真工具，ANSYS一直讲多物理场仿真作为支撑和核心技术。尤其以Workbench作为多物理场仿真的平台。每完成一次仿真软件收购之后，就会将其融入到多物理场平台中。

功能较全能提供多物理场耦合仿真的软件有ADINA，Altair Multiphysics，

STAR-CCM ，其它号称能解决多物理场问题的工具其实都算不上真正意义上的多物理场工具，比如MSC的SimXpert，达索的SIMULA等。

后记

在仿真软件（CAD/CAE/EDA/CFD）开发工作中，发现很多软件研发测试人员没有相关背景 希望该篇文章能有所帮助。