1.1 PDE
# y& ]8 C R- a, V2 T有限元是一种求解问题的数值方法，求解什么问题呢？－－求解PDE(偏微分方程).那么PDE是做什么用的呢？－－描述客观物理世界。我想如果这两个问题搞清楚了也就明白了为什么要用fem,fem可以做那些东西。 PDE可以描述很多物理现象，电磁，流体，换热，diffusion,力学，河床变迁，物种竞争，股票金融，等等等等。。。。乃至整个宇宙，当然也不是所有的物理现象都可以用PDE描述，如微观世界分子原子的运动等等，所以我从来都不建议用有限元方法仿真微观物质现象的原因，这是题外话就不多说了。除了PDE以外，ODE同样也可以描述客观世界，但ODE多用于控制系统，很多线性PDE的解法也都是将PDE转化为ODE来解的，如separating variables方法。
1.2 求解PDE
5 g. U7 J' ?# y2 ]$ a\! |/ Y5 ~/ n( V( V) p0 }' ^
有了PDE以后，问题是如何求解并得到结果，首先要说明的是不是所有的PDE都有解的，往往有解的PDE才有实际工程意义。对于数值解法，常用的是有限差分，有限元和谱方法，还有蒙特卡罗法。有限差分出现的较早，计算精度相对较高，但是费时，且模型形状必须规则，边界条件处理困难。有限元方法效率高且满足精度要求，边界条件容易处理，得到了广大的应用。谱方法是基于FFT方法的解法，精度高，收敛快，必须使用periodic boundary condition，适合微观尺度的PDE解。蒙特卡罗法不是基于弱解形式的，多用于金融分析，这里还是着重有限元解PDE，顾名思义，有限元将整个计算几何模型划分为很多小的单元（element）,每个单元的含有一定数量的节点(node)，具体单个单元有多少节点，有对应的不同算法与差值方程，拿一个简单的线性4节点平面单元来说，每个单元包含4个节点，每个节点有对应的variable值，比如简单固体力学问题，每个节点就有对应的位移值，热力学问题每个节点就有对应的温度值，等等。然后单元内部的variables就通过差值方法计算得出。
* B( P+ k6 N+ e% B
I# \\# ^- p. P {1.3 Galerkin approximation 与 weak formation
/ t' U0 w# _* v4 G/ A* w- M$ @3 w, F
弱解(weak formation)是建立在变分法基础上的，通过这个方法将strong form的PDE转换为weak form,使得有限元的求解成为可能，具体如何推导weak form可以参考一些有限元书籍，如果一本基础的有限元书籍没有介绍如何推导weak form,那么可以考虑选择换一本了。推导所得的弱解形式仍需要通过计算机来计算，Galerkin方法推导出含有求和符号的公式，在计算机中多半以loop的形式来计算这个量，往往这个量就是stiffness matrix中的component。公式中还会存在积分计算，有限元方法多用gauss quaradure的方法来计算，精度一般可以满足。也就说一般有限元计算中存在两次approximation,一次是Galerkin一次是gauss,这也是很多人在计算完以后需要进行validation的原因。单个单元的stiffness matrix计算完成后，还需要将所有单元的矩阵装配为一个大型的矩阵，然后进行线性代数计
算。这个装配是很有技巧的，因为一般情况下stiffness matrix是一个很大的稀疏矩阵，０值往往可以省去以节省计算量。要知道，一个１００个单元的４节点平面单元单个变量的stiffness matrix会有101x101大，随着单元数或变量的增加，计算是惊人的。
$ b0 Z- w3 N\3 p3 p7 x
1.4 后处理
其实对于最基本的有限元方法，求解得到的仅仅是variable的值，如力学就是节点位移值，thermal problem就是温度值，流体就是位移速度加压力值，如果我们想知道应力或者应变怎么办呢？后处理系统里面个都会增加相应子程序来计算stress, strain, flux等等。这也就是为何我们能看到各种各样俄contour的原因了，当然读者也可以自己加入计算各种量的子程序，如应变能密度什么的。关于什么是有限元就介绍到这里，仅仅是一些随笔和想法，具体的理论和推导需要自身实践与探索，本文行文仓促只是阐述自己对有限元的粗浅理解。有不对的地方还请指正。下一章会谈及一些我曾经用过得软件。 CH2_有限元软件的介绍与比较 by caoer
\
有限元软件很多，有商业的，开源的，免费的，并行的，多物理场的，专业于某领域的，这里仅仅介绍一些笔者曾经用过的，遗漏的还请高手补充，我会注明是你的补充。
+ M. }% Z( l' H, [; ]% m1 D- T. T' P5 U& O
7 ?6 N; Q- _1 S2.1 ANSYS
第一个学习的有限元软件，也是上世纪末本世纪初国内最流行的，为什么流行呢？应为ansys的教程铺天盖地，所以大家都学习ansys因为有教程上手快，那个时abaqus的教程可谓是凤毛麟角，自然学习的人也就不多，当时组里导师都是用ansys，自然我也就用了。如果一个有限元软件推广的好，那么他的教程一定要推广的好，fem是一个专业性很强的软件，教程推广不够，销售自然不行。题外话，接下来说说对用他的感受。
4 @& b5 N' h G: d) T# b0 Z) `
当时学习ansys是直奔apdl去的，加上理论背景很弱所以学习过程异常艰苦。总得说来ansys是一个很中规中具的软件，功能很强大，计算也很可靠，速度快，基本上固体力学的问题都能解决。作为一个工程软件还是很不错的。帮助文件很丰富，我也挑不出他什么缺点。只是感觉国内学习ansys已经陷入了一个怪圈，用户喜欢就某一个应用技巧钻的很深，对于其基础的理论却忽略不视，这无不于大量ansys教材的编写误导有关。是指算例，却不关心原理，悲哀。。。
2.2 ＣＯＭＳＯＬ （femlab）
: _& K) }& F; [/ U( O7 B \\+ D0 Y' g6 ]1 c\
第一次用comsol的时候它叫femlab2.0,直奔多物理场耦合而去的，无奈的是只是作为一种兴趣学习的，因为时间精力问题最后也没有很深入，偶然机会１年后又开始认真的学习，逐渐体会到他的强大，我想说的是comsol对于各种物理场的research的确是一把利器，他可以任意的由用户输入ＰＤＥ,计算结果并出contour,这对于懂得有限元理论却又不想花大量的时间来编程的人来说，就是种bliss,然而缺点就是速度太慢，这也是可以理解的，毕竟
comsol是将单元与pde信息进行重组后计算，而不是很多商业有限元已经将Ｐde固化在单元信息里了，但是其出色的weak form PDE模块是不可多得的有限元分析模块。还有其建模与后处理模块相当方便，比ansys方便很多。至于计算的精确性笔者还抱有怀疑态度，细节就不提了。国内现在这个软件推的很好，大有赶超ansys之势。
% V8 \\& _5 v) T' d+ `7 `# m ?7 N
9 ?0 E# Q+ Q& B, P' e2.3 FEAP
这个其实是老大哥，虽然很多人不知道它。开源，不免费，但是很便宜，也有免费的学生版本，大家可以下载学习用。据说abaqus和ansys都是源于这套程序，abaqus技术部里面的人都是feap搞这个组里出来的，虽然是用fortran编的，但是编程思路清晰，注解丰富，参考文档相近，实在是不可多得学习fem的神器。不提供gui的输入方式，命令的输入和apdl和umat方式很像，都是文本macro方式。计算速度超快，带有后处理功能，提供强劲的二次开发接口，可以自己编写子程序输出其他后处理软件如tecplot.具有openmpi并行能力，提供丰富的elastic, plastic，hyperelastic 模块。缺点是之用于固体力学与热计算，当然可以自己开发其他单元，如电磁的。还有就是feap只运行于linux下，使用者需要知道如何make程序。
% \\( {+ m5 V' d, {. m
2.4 libmesh 和 deal II
之所以把这两个程序放在一起，是因为都出自于同一个大师。c++开源免费程序，非常的robust，但是还有一些小bug存在，可以应用于更广泛的领域如流体，强劲在于adaptive meshing,和remesh,从而可以解决很多常规有限元无法解决的问题，如singularity问题。需要有c++背景以及有限元理论知识，也是笔者遇到门槛最高的一个有限元软件，推荐给专业人事。同样是在Ｌinux下运行的软件。
2.5 ALGOR_pipeline
h4 Q1 n. B8 z) b1 N这个软件在也挺有名，但是国内用的少，主要是固流耦合这快做得比较好，就我所知一些有名的飞机和流体机械厂都用他，我主要是用pipeline这个插件，也许名字记得不对，就是用来计算整个管路系统的震动的，前处理中设定管路，就可以计算出各阶振型与频率，很方便，据说准确性挺好。其主模块用的不多，不做讨论了。
基本上笔者主要使用的就是这几种有限元软件，我想如果能精通上面任何一种软件，学习其他的软件都是一个很快的过程。如果你是博士生且有一个比较长得研究学习时间，推荐用feap和deal II,如果是硕士，推荐comsol，本科毕业设计，推荐ansys.当然每个人水平周围环境都不一样，不能一概而论，总之多学一点东西没有坏处，每个软件都有其深厚的背景。也欢迎高手补充其他软件，投条给我。 yCH3_有限元的数值方法 by caoer
( S. x! y9 C5 h+ [
这里先略过有限元的几何建模与网格划分部分，因为一来不是数值方法的主要部分，二来我对这块也不是很精通，所以直接从数值解法入手。这部分可以说是ch4与ch5的总体概览部分。也是fem的核心。
4 U. W3 b% ?' a( s! t6 i
3.1 单元离散化与jacobian
划分好网格后，就意味着单元编号以及节点都已确定下来了，其实这步在有限元分析的一开始就设定好了。拿平面单元为例，如果选用平面8节点矩形单元来计算热力学问题，那么每个单元有8个节点温度值，要知道每个节点的位置并不是规则均匀分布在实际有限单元上的，人们通常将节点的global坐标转换到local坐标上，以方便计算推导，之后再通过jacobian转换到global上，这点和连续介质力学的reference转换是同一道理。
6 m) X& p R$ {
/ [) |/ E( _4 v+ X. J% D3.3 运动方程与各种矩阵
得到了单个单元的8节点test function 或 grediant of test function的积分值以后（积分的计算用gauss法），就需要联系所有的单元，装配成为一个整体的矩阵，这个矩阵就是stiffness matrix,如果是一个4单元简单正方形区域，那么装配好的stiffness matrix就是一个21x21矩阵，以k标记，除此以外，如果是瞬态问题，也会有质量矩阵存在，多半是一个对角矩阵，其值一般是test function的积分值，如果存在damping那么还会存在damping matrix,一个有限元系统一般只存在这3个矩阵，然后通过一个运动方程将其连立起来，即
M*u_tt+D*u_t+K*u=F,F是source。这就是将pde转化为弱解再转化为有限元方程的最后形式。
* x) X7 [6 }) @+ U( v2 l
3.4 Ax=b求解
有了 M*u_tt+D*u_t+K*u=F 方程，接下来要做得就是求解他，如果是一个2x2的矩阵，手工即可计算得到，实际应用中，往往都是上万阶的矩阵。也许有人会问了u_tt和u_t是怎么解得的？瞬态问题由于有了时间变量的加入，所有必须要有对应的时间积分求解器，这点在求解器部分会详述，经过时间求解器的计算后，运动方程还可以化简为 K_tilde＊u=F的形式，也就是线代求解器需要解的最后方程，求解后，各节点的u值得到。完毕。
6 ]9 |3 H' M7 Y4 }
3 j# `8 K& X, F8 J0 o% W6 3.5 多场耦合方法
如果有多个场存在怎么办呢？道理很简单，比如一个固体场一个温度场，两个pde出来的运动方程是这样的 M1*u_tt+D1*u_t+K1*u=F1和 M2*u_t+K2*u=F2,将两个场的运动方程线性进行叠加（线性耦合），得到 M1*u_tt+(D1+M2)*u_t+(K1+K2)*u=F1+F2 ,再化简，得到 K_tilde＊u=F3，使用求解器求解既得结果。这里只是简单的线性耦合，对于复杂物理现象的非线性耦合都是在这个基础之上进行变化。而且要注意的是，最后所得的K_tilde不一定是对称矩阵，也就意味着求解器必须要有解Unsymmetric matrix的能力。
7 v! k5 E- y' K% p! R1 F$ K) G\
3.6 边界条件加载
边界条件的加载，往往都是在单元矩阵装配完毕以后进行的，往往是variable或是gradiant of variable的值被限定住了，如果pde存在3阶grediant，那么还会存在laplace operator of veriable的边界条件。具体编程的实现我能完全总结出来，先欠着回头再写，以免误人子弟。
1 f. F: }5 D% G# m* e0 ~
帖子是想到什么写什么，难免凌乱没有章法，也有不少遗漏之处，想到了再回来补 CH4_常用固体单元 by caoer
$ N! V4 K! X; y* U: y7 s/ f, D9 @\
记得自己刚学ansys的时候就被自带的100多个单元模型搞的晕头转向，到底应该选择什么样的单元？这些单元有什么不同？这些问题一直困扰着我，其实把每一种单元都搞透是一件不太可能的事情，但有件事情可以做得就是对于几个大类的单元有个底层数值方法上的了解，这对于单元的选择与开发是至关重要的。这里仅以3d单元作为例子。
; L% \\6 s. o) X
有限元方法中以单个单元为基础，所有的单元都有相同的属性，这些单元组成了整个domain,可以说了解了单元也就了解了整个模型，现代有限元软件都将pde整合在了单元信息里面，还句话说，Pde所推导出的弱解形式中非边界条件term都体现在了单元中。
4.1 固体单元, 变量：位移u1,u2,u3
+ d3 W; I: X* n' W工程力学中最常用的单元，常用于连续介质的力学计算，一般分为small deformation 和 finite deformation 两个区域，记得ansys里面有个大变形开关，一打开后计算变得惨不忍睹，所以选择small还是finite要看具体的问题，常见工程材料一般都是small deformation,此外对于不同的材料，如超弹，粘弹，弹塑材料都有不同的数学模型，这里选择的时候一定要慎重，不要想当然，最好要看看手册和理论，计算方法上常用的有
displacement,mixed,enhanced strain方法，displacement是最原始的，收敛效果不好，可以用于解决一般基础问题，推荐后面两种来解非常规材料。此外，对于同样的pde模型，单个单元节点越多计算越精确，单元数目相同情况下，八面体单元比6面体单元精度高。
4.2 热单元，变量：温度T.
是相对简单的一种单元，pde本身也很简单，就是换热方程，节点自由度就是T,由于温度T没有方向性，所以T只有1维量，计算量小。另外heat flux = gradient of T, 重点参看一下Fourier heat conduction equation.
4.3 壳单元，变量：u1 , u2 , u3 , θ1 , θ2 , θ3
有些具有曲面外壳的材料，其一个方向的变形要远远小于其他方向的，如人的骨头表面，就是一层非常硬的皮质骨，所以生物力学的人都用shell单元作表面单元.也有一种大变形的壳单元，只含有5个变量,少一个θ或是u,但是有可能在计算的时候不收敛，所以在使用不同的shell单元时，一定要多看看说明，遇到不收敛的情况也不用慌张，多从手册分析解决问题。