本文摘要:(由ai生成)
本文总结了作者自2003年以来学习有限元方法的经验。介绍了有限元法的基本概念、多款有限元软件的使用体验、数值方法深入解析及固体单元应用。探讨了直接、迭代、时间和Newton求解器。展望了有限元法的未来,强调并行计算和多尺度模型计算的重要性。作者感谢有限元社区并鼓励大家为领域贡献。
接触有限元时间不算常也不短,从03年到现在(2010年)吧,中间断断续续的学习有限元,大大小小软件也用了不少,牛人的帖子大师的著作读了一些,发现还有很多学弟学妹初学fem,提出的问题五花八门但大多不对正路,不禁回想起了当年自己学习有限元时情景,孤立无助到处瞎碰,通宵熬夜啃手册也得不出一个所以然来,现在回头想想如果有几个良师益友会是多么幸福的事情,写此文仅仅只为了让后生少走些当年自己走过的弯路,有更多的精力去研究更深层次的问题,而不是在门槛上打转。
本文并不是一套完整的教科式学习材料,对于真正想在有限元上做一些工作的,还是建议多读大师的著作和程序,少看此文,另外,此文仅仅是自己对有限元法粗浅的理解,如有疏漏的地方,还请各位大师补充指正。
1. PDE
有限元是一种求解问题的数值方法,求解什么问题呢?求解PDE(偏微分方程)。那么,PDE是做什么用的呢?描述客观物理世界。我想如果这两个问题搞清楚了也就明白了为什么要用有限元,有限元可以做哪些东西。PDE可以描述很多物理现象,电磁、流体、换热、声学、扩散、相变、各种力学、河床变迁、物种竞争、股票金融等等等等......乃至整个宇宙。当然,也不是所有的物理现象都可以用PDE,所以我一般不建议用有限元方法仿真微观物质现象的原因,但也有PDE应用于微观物质并得到很好的结果,如泊松方程来解析plasma的物理现象,这在量子物理里用统计的方法过于庞大,泊松方程反而使问题简单而且能吻合实验,这些都是题外话就不多说了。除了PDE以外, ODE同样也可以描述客观世界, ODE相对简单,多用于控制系统,有一些线性PDE的解法也都是将PDE转化为ODE来求解的。
2. 求解PDE
有了PDE以后,问题是如何求解并得到结果。首先要说明的是,不是所有的PDE都有解的,往往有解的PDE才有实际工程意义。对于数值解法,常用的是有限差分、有限元、有限体积和谱方法,还有蒙特卡罗法。
有限差分出现地较早,计算精度相对较高且可控,但模型形状必须规则,边界条件处理困难;好处是可以比较方便的控制计算精度,编程简便,固定节点的网格划分形式适用于流体类的仿真。
有限元方法效率高且满足精度要求,边界条件容易处理,得到了广大的应用,尤其是在固体领域。
有限体积类似是有限差分和有限元的结合,特点是计算点在单元中心,雷曼边界处理单元之间的交接,处理流体中的激波有着特殊优势,计算方法边界、速度要比同等级的有限元方法快。
谱方法由于可以采用高阶差值方程和FFT方法来求解,使得程序有着精度高,收敛快的特点,也克服了有限元条件下使用高阶插值方程计算费时的缺点,常常使用periodicboundarycondition,但也有越来越多的算法使得一类二类边界成为可能,适合微观尺度的PDE解,谱方法和有限元结合产生的谱元法取两者之优点,使得应用前景非常好。
蒙特卡罗法不是基于弱解形式的,随机数的多维采样最终得到统计上的结果,多用于金融分析。
咱这里还是着重有限元解PDE,顾名思义,有限元将整个计算几何模型划分为很多小的单元 (element),每个单元含有一定数量的节点 (node),具体单个单元有多少节点,有对应的不同算法与差值方程。拿一个简单的线性4节点平面单元来说,每个单元包含4个节点,每个节点有对应的自变量值,比如简单固体力学问题(位移法),每个节点就有对应的位移值,热力学问题每个节点就有对应的温度值,等等。然后,单元内部的结果就通过差值方法显示出来。
3. Galerkin approximation与 weak formation(变分/虚功)
弱解 (weakformation) 是建立在变分法基础上的,通过这个方法将strong form 的PDE转换为weak form,使得有限元的求解成为可能,具体如何推导weak form可以参考一些有限元书籍,如果一本基础的有限元书籍没有介绍如何推导weak form,那么可以考虑选择换一本了。
推导所得的弱解形式仍需要通过计算机来计算, Galerkin方法推导出含有求和符号的公式,在计算机中多半以loop的形式来计算这个量,往往这个量就是stiffness matrix 中的各个量。公式中还会存在积分计算,有限元方法多用gauss quaradure的方法来计算,精度一般可以满足。也就说一般简单的限元计算中存在两次approximation,一次是Galerkin,一次是gauss 法求积分,这也是很多人在计算完以后需要进行验证的原因,同时对于非线性问题newton 法本身就是通过计算误差来判断收敛的,固体有限元常常使用能量增量活residual 增量作为newton 求解器的判据。单个单元的stiffness matrix计算完成后,还需要将所有单元的矩阵装配为一个大型的矩阵,然后进行线性代数计算。这个装配是很有技巧的,因为一般情况下stiffness matrix 是一个很大的稀疏矩阵,0值往往可以省去,以节省计算量。一个由10×10×10个8节点矩形单元组成的模型会有11×11×11个 node,如果是热力学问题自变量是温度T,每个节点上有一个自由度,组装后的的stiffness matrix 会有 (11×11×11)×(11×11×11)=1771561之大。
随着单元数或变量的增加,计算是惊人的,这样一个大的矩阵求解显然不能用常规的方法(如gaussian elimination法),这就是各种这样求解稀疏矩阵算法存在的原因了,稀疏矩阵存储转换方法可以参见Bathe的《FiniteElementProcedure》,有效且快速求解线代方程组是一个好的solver 的标准
4. 后处理
其实对于最基本的有限元方法,求解得到的仅仅是各节点的自变量的值,对于基本的位移法,如力学就是节点位移值,热力问题就是温度值,流体就是位移速度加压力值。如果我们想知道应力或者应变怎么办呢?后处理系统里面个都会增加相应子程序来计算stress,strain,flux等等。这也就是为何我们能看到各种各样的云图的原因了,当然读者也可以自己加入计算各种量的子程序,如应变能密度什么的。
关于什么是有限元就介绍到这里,仅仅是一些随笔和想法,具体的理论和推导需要自身实践与探索,本文行文仓促只是阐述自己对有限元的粗浅理解。有不对的地方还请指正。接下来会谈及一些我曾经用过的软件。
有限元软件很多,有商业的、开源的、免费的、并行的、多物理场的、专业于某领域的,这里仅仅介绍一些笔者曾经用过的,遗漏的还请高手补充。
1. ANSYS
第一个学习的有限元软件,也是上世纪末本世纪初国内最流行的,为什么流行呢?因为ansys的教程铺天盖地,所以大家都学习ansys因为有教程上手快,那个时abaqus的教程可谓是凤毛麟角,自然学习的人也就不多,当时组里导师都是用 ansys,自然我也就用了。如果一个有限元软件推广的好,那么他的教程一定要推广的好,有限元是一个专业性很强的软件,教程推广不够,销售自然不行。题外话,接下来说说对用他的感受。
当时学习ansys是直奔apdl去的,加上理论背景很弱所以学习过程异常艰苦。总的说来ansys是一个很中规中具的软件,功能很强大,计算也很可靠,速度快,基本上固体力学的问题都能解决。作为一个工程软件还是很不错的,帮助文件很丰富,我也挑不出他什么缺点。说到Ansys,不得不提到他的创始人Swanson博士,据说物理出身的他老人家当年每时每刻都在写程序推公式,飞机上、睡觉前,永远都在工作,这种对有限元的坚持永远是我们学习的榜样。说到Swanson不得不提到Swansea,后者由于大Z的存在,可以说是有限元方法的发源地,可惜现在这个Swansea大学逐渐没落了,大Z也在09年的时候去世了,留下了《 Finite Element Method》第六版的绝唱。回到ansys,现在这个公司的技术核心似乎是老印当家, 这也是我对其将来发展不太看好的原因之一。
2. COMSOL
第一次用comsol的时候它叫femlab2.0,无奈的是只是作为一种兴趣学习的,因为时间精力问题最后也没有很深入,偶然的机会在1年后又开始认真的学习,逐渐体会到他的强大,我想说的是comsol对于各种物理场的研究的确是一把利器,他可以任意的由用户输入PDE,计算结果并出云图,这对于懂得有限元理论却又不想花大量的时间来编程的人来说,就是种天赐之福。建模和界面设计非常人性化,上手便捷,然而缺点就是速度太慢,这也是可以理解的,毕竟comsol是将单元与PDE信息进行重组后计算,而不是像很多商业有限元已经将PDE固化在单元信息里了,但是其出色的weak form PDE模块是不可多得的有限元分析模块。还有其建模与后处理模块相当方便,比ansys方便很多。至于计算的精确性笔者还抱有怀疑态度,细节就不提了。
国内现在这个软件推的很好,大有赶超ansys之势。值得一提的是comsol的并行计算能力,采用的umfpack+mpich的构架,总体并行效率还是非常令人满意的。由 于umfpack到5.2版本后就改变了软件协议,现在的comsol应该不再采用最新的 umfpack了,comsol在招聘的时候也只是要求blas,c++和mpi,可见comsol是一个会不断发展的软件。
3. FEAP
这个其实是老大哥,虽然很多人不知道它。开源,不免费,但是很便宜,也有免费的学生版本,大家可以下载学习用。据说abaqus和ansys都是源于这套程序,abaqus技术部里面的人都是feap搞这个组里出来的,虽然是用fortran编的,但是编程思路清晰,注解丰富,参考文档相近,实在是不可多得学习有限元的神器。不提供gui的输入方式,命令的输入和apdl方式很像,都是文本macro方式。计算速度超快,带有后处理功能,提供强劲的二次开发接口,可以自己编写子程序输出其他后处理软件,如tecplot。具有mpi并行能力,提供丰富的elastic、plastic、hyperelastic模块。缺点是只用于固体力学与热计算,当然可以自己开发其他单元,如电磁。虽然其并行化在8.1版本的时候就已经实现,但基于PETSc的并行求解方式改变了FEAP独立运行的最大特色,大多数用FEAP的人应该并不打算做大规模计算的。
FEAP俨然已经成为了大量有限元的开发蓝本,我就看到过有好几款开源的有限元软件和FEAP有着极为相似的编程手法与思路。另外,就是feap只运行于linux下,使用者需要知道如何make程序。
4. libmesh和deal II
之所以把这两个程序放在一起,是因为都出自于同一个大师,一位deal II的主创人员从德国来到了德州,组织开发并行且网格自适应的libmesh。c++开源免费程序,非常的robust,可以应用于更广泛的领域,如流体。强劲在于adaptive meshing和并行,从而可以解决很多常规有限元无法解决的问题,如 singularity问题,大规模问题。libmesh虽然是后起之秀,但是发展异常迅猛,无奈只有那20多个例子作为教程,对于初学有限元的朋友那就是灾难了,而且似乎固体不是libmesh的主攻方向,如果要在libmesh上开发一个壳单元,就必须对整 个内部插值环境有很好的了解,只能说难度很大。需要有c++背景以及有限元理论知识,是笔者遇到门槛最高的一个有限元软件,推荐给专业人事。同样是在Linux下运行的软件。
Libmesh最大的优势应该是大规模问题,应用对象也是cluster。
5. ALGOR_pipeline
这个软件在也挺有名,但是国内用的少,主要是固流耦合这快做得比较好。据我所知,一些有名的飞机和流体机械厂都用他,我主要是用pipeline这个插件,也许名字记得不对,就是用来计算整个管路系统的振动的,前处理中设定管路就可以计算出各阶振型与频率,很方便,据说准确性挺好。其主模块用的不多,不做讨论了。
基本上笔者主要使用的就是这几种有限元软件,我想如果能精通上面任何一种软件, 学习其他的软件都是一个很快的过程。如果你是博士生且有一个比较长得研究学习时间,推荐用feap和 libmesh或者是基于如PETSc/blas求解器这样的自编程序。如果是硕士,推荐omsol。本科毕业设计,推荐ansys/abaqus/adina等。
这里先略过有限元的几何建模与网格划分部分,因为一来不是数值方法的主要部分,二来我对这块也不是很精通,所以直接从数值解法入手。这部分可以说是第3章和第4章的总体概览部分。也是有限元的核心。
1. 单元离散化与Jacobian
划分好网格后,就意味着单元编号以及节点都已确定下来了,其实这步在有限元分析的一开始就设定好了。拿平面单元为例,如果选用平面8节点矩形单元来计算热力学问题,那么每个单元有8个节点温度值,要知道每个节点的位置并不是规则均匀分布在实际有限单元上的,人们通常将节点的global坐标转换到local坐标上,以方便计算推导,之后再通过jacobian转换到global上,这点和连续介质力学的 reference转换是同一道理。
但是,jacobian的建立是不是必须的呢?不是,如果所有的插值方程都是基于 global的,那么local空间可以省去,jacobian自然不需要了,但这种做法不被推荐,应该编出来的程序可扩展性非常差,而且计算量较大,所以在大型通用有限元里面都是计算jacobian的。Jacobian做为构建两个坐标系关系的偏微分矩阵,其值必须要大于 0。
2. 运动方程与各种矩阵
得到了单个单元的8节点形函数和形函数导数的积分值以后(积分的计算用gauss法),就需要联系所有的单元,装配成为一个整体的矩阵,这个矩阵就是stiffness matrix,如果是一个4单元简单正方形区域,那么装配好的stiffness matrix就是一个 21×21矩阵,以k标记。除此以外,如果是瞬态问题,也会有质量矩阵存在,多半是一个对角矩阵,其值一般是形函数的积分值,如果存在damping,那么还会存在damping matrix。一个有限元系统一般只存在这3个矩阵,然后通过一个运动方程将其连立起来,即M*u_tt+D*u_t+K*u=F,是source。这就是将pde转化为弱解再转化为有限元方程的最后形式。这个矩阵绝对是一个稀疏矩阵,对于稀疏矩阵传统的求解方法似乎就显得非常笨拙了(如blas、MKL),所以各种各样求解稀疏矩阵的求解器就出现了 (gotoblas,PETSc)。
3. Ax=b求解
有了M*u_tt+D*u_t+K*u=F方程,接下来要做的就是求解他,如果是一个2×2的矩阵,手工即可计算得到,实际应用中,往往都是上万阶的矩阵。也许有人会问u_tt和u_t是怎么解得的?瞬态问题由于有了时间变量的加入,所有必须要有对应的时间积分求解器,二阶的有newmarks、HHT、energy conserved方法,一阶的有crank nicolson,向前向后Euler法等等。这点在求解器部分会详述,经过时间求解器的计算后,运动方程还可以化简为K_tilde*u=F的形式,也就是线代求解器需要解的最后方程,求解后,各节点的u得到。
4. 多场耦合方法
如果有多个场存在怎么办呢?道理很简单,比如一个固体场一个温度场,两个pde出来的运动方程是这样的M1*u1_tt+D1*u1_t+K1*u1=F1和 M2*u2_t+K2*u2=F2,将两个场的运动方程线性进行叠加(线性耦合),得到 M1*u3_tt+(D1+M2)*u3_t+(K1+K2)*u3=F1+F2,再化简,得到 K_tilde*u3=F3,使用求解器求解既得结果。这里只是简单的线性强耦合,对于复杂物理现象的非线性耦合都是在这个基础之上进行变化。而且要注意的是,最后所得的K_tilde不一定是对称矩阵,也就意味着求解器必须要有解Unsymmetric matrix的能力。此外,还有各种各样的弱耦合方法,如果是自编程序的话,弱耦合是非常多样的,这取决于程序设计的意图和目的。
5. 边界条件加载
边界条件的加载,都是在所有单元矩阵装配完毕以后进行的,往往是自变量或是自变量梯度的值被限定住了,对于低阶边界条件,如 Dirichlet boundary condition,常用的方法是增加一个惩罚量,这个惩罚量相当大>1e10,这样便使得u被固定在人为设定的u0位置。对于natrual boundary condition,也就是q_n的值如何确定,通常与变量的梯度值相关,如果q_n的表达式已知,那么只需要计算出q_n的精确值带入source即可。需要注意的是,natural boundary condition要求比较精细的网格,否则会导致局部节点的计算失真。对于固体问题,q_n代表的是压力,具体的节点力(或弯矩在板壳杆单元中)可以通过修改源相{F}来施加载荷。
记得自己刚学ansys的时候就被自带的100多个单元模型搞的晕头转向,到底应该选择什么样的单元?这些单元有什么不同?这些问题一直困扰着我,其实把每一种单元都搞透是一件不太可能的事情,但有件事情可以做,就是对于几个大类的单元有个底层数值方法上的了解,这对于单元的选择与开发是至关重要的。这里仅以3d单元作为例子。
有限元方法中以单个单元为基础,所有的单元都有相同的属性(也可组合各种单元),这些单元组成了整个区域,可以说了解了单元也就了解了整个模型,现代有限元软件都将pde整合在了单元信息里面,换句话说,Pde所推导出的弱解形式中非边界条件的元素都体现在了单元中。
1. 固体单元,自变量:位移u1、u2、u3
工程力学中最常用的单元,常用于连续介质的力学计算,一般分为small deformation和finite deformation两个区域,记得ansys里面有个大变形开关,一打开后计算变得惨不忍睹,所以选择small还是finite要看具体的问题。此外对于不同的材料,如超弹、粘弹、弹塑材料都有不同的数学模型,这里选择的时候一定要慎重,不要想当然,最好要看看手册和理论,计算方法上常用的有displacement、mixed、enhanced strain方法,displacement是最原始的,收敛效果不好,可以用于解决一般基础问题,推荐后面两种来解非常规材料。此外,对于同样的pde模型,单个单元节点越多计算越精确,单元数目相同的情况下,6面体单元比4面体单元精度高。
商业软件一般都考虑的很周全,如果是自编软件的话,就需要考虑很多问题了,如locking、buckling、convergence等等,这些东西很有趣也很深奥,既需要深厚的功底,又要有实际的编程能力,如果要做实际可用的大规模问题,对计算机CPU的构架、cache的大小、data layout和编译器都需要很好的理解,才能编出高质量的代码。所以,编写一个small deformation的固体单元对于科技发展的今天,已经是本科生的课堂设计实验课了。但对于大变形、多场耦合、超弹/粘塑性材料、多尺度,还是需要比较深厚的张量推导内力的,所以搞力学理论背景的要多了解计算机相关原理及建立高效编程习惯,搞计算机出身的要建立张量运算推导和物理现象的概念,那样就是有限元开发的有用人才,这句话我和广大有限元爱好者共勉。
2. 热单元,变量:温度T
是相对简单的一种单元,pde本身也很简单,就是换热方程,节点自由度就是T,由于温度T没有方向性,所以T只有1维量,计算量小。另外,heat flux=gradient of T,重点参看一下Fourier heat conduction equation方程。
3. 壳单元,变量:u1、u2、u3、θ1、θ2、θ3
有些具有曲面外壳的材料,其一个方向的变形要远远小于其他方向的,也有一种大变形的壳单元,只含有5个变量,少一个θ或是u,但是有可能在计算的时候不收敛,所以在使用不同的shell单元时一定要多看看说明,遇到不收敛的情况也不用慌张,多从手册分析解决问题。
壳单元很多,最著名的要属MITC了,Bathe在84年开发的,影响至今还在,基本上目前所有的壳单元(除了degenerated solid method)都是基于MITC的。壳单元(理论)可以说是近100年来固体力学工作者最令人神往的地方,因为其广泛的应用使得大家不得不对其感兴趣,小到人的眼球,大到屋顶,重要到导弹潜艇,到处都是壳啊……壳太复杂了,以至于有很多人尽量避开他去研究板,要说目前最高深的板壳理论还是钱伟长大师的非线性板壳理论(个人理解),可惜在现在西方的文献里,大家都只提von karman nonlinear plate/shell theory,钱老的研究可以说是提前的人类文明至少50年(还是个人理解),因为就目前来说nonlinear板壳的应用还是凤毛麟角。将nonlinear shell/plate理论变成有限元单元的,目前笔者还没看到,等待某位大牛出现将其编出。
4. 框架单元,变量:u1、u2、u3、θ1、θ2、θ3
常用于结构计算,如弯曲、剪切变形。也就是材料力学里面各种梁的计算,这里pde方程是关于u的4阶梯度,解的方法是将自变量的二阶梯度看作是Moment,3阶看作是shear stress,这类单元通常是2节点单元,一头一尾相互连接。
索单元和板单元就不赘述了,可以参看相关资料,与框架单元大同小异。索单元甚至更简单,和热单元类似。还有一些薄膜单元点单元什么的应用很少,就不多说了。
看到simwe上很多人开口闭口二次开发,其实所谓二次开发很大程度上都集中于新单元的开发,如果你有新的pde或者是本构,要想计算实现它,就必须写出新的单元代码,从而实现计算。其他方面的二次开发无外乎一些后处理或是mesh上面的新功能,而这些商用有限元基本成熟,所以再提二次开发之前,最好先考虑现有的单元库里面有没有我需要的单元,如果没有再二次开发它,而这个过程是漫长而艰苦的。
最后的问题都围绕求解Ax=b这个方程组了,也就是要求解它,那么高效的求解器是必要的。其实很多商业有限元的使用者都忽略求解器的选择,实际上求解器直接决定了求解的精确性和速度,选择正确有效的求解器也是所有分析人员必须了解的问题。有人抱怨有限元计算太费时或是老报错,很有可能选择了不正确的求解器,这里仅对主流一些求解器作一些介绍。毕竟solver成千上万,很多人也喜欢编写自己的求解器来解决特定的问题,这样他们认为更加有效率且误差小。本节主要介绍线代求解器和时间求解器,线代求解器又分为直接求解和迭代求解器。
1. 常用直接求解器
最经典的且唯一的直接求解器非gaussian elimination莫属了,pivoting的引入使得这种方法的稳定性大大增强。虽然此方法易于并行,但对于大规模问题,其计算速度是成倍的下降而内存的使用是成倍的增加,所以现代大型有限元基本放弃了它。直接法中Factorization是一个必要的过程,Factorization的快慢好坏也决定了求解成果。LU Factorization有着很好的计算效率和通用性。对于对称正定矩阵A,Cholesky Factorization是一个高效的分解方法,记得前阵子simwe还有朋友问Cholesky如何发音,可见这个方法还是很流行的。QR Factorization对于求解特征向量系统是非常有效的,但相对开销较大。我认为直接求解器仅仅是指gaussian elimination方法而已,只不过在Factorization和backsubstittution上有会有些不同的分支,其总原理是基本一样的。
2. 常用迭代求解器
迭代求解器在求解效率上有着惊人的提升,加上各种改良的算法,如relaxation,incomplete factorization,multigrid等一系列加速方法,使得应用非常广泛。但迭代求解方法可能会有求解器不收敛的情况发生,这是在求解具体问题时需要避免的。基本的迭代法应该是Jacobi方法了,编程非常容易,Gauss-Seidel是一 种类似的方法,relaxation的引入使得相对小的迭代次数就可以得到的收敛的结果。Jacob和Gauss-Seidel都是非常易于并行化的迭代方法。对于非正定或非对称刚度矩阵,还有一些细节处理使得系统半正定。Preconditioner在迭代法中有着重要的作用,本质就是在线代方程两边乘以一个矩阵的逆,P^{-1}Ax=P^{-1}b,其中矩阵P就是preconditioner,选择一个好的矩阵P可以使系统的特征值远离0。如何选择好的P矩阵,有很多相关书籍和论文介绍大家可以自己去看。Preconditioned Conjugate Gradients(PCGs)方法由于利用了之前的计算结果,可以大幅提高计算速度,这种方法牺牲了存储提升了速度,是基于Krylov子空间的经典方法。还有一些基于Krylov空间的求解方法,如大名鼎鼎的General minimal residual(GMRES),其本质是multigrid和Krylov子空间方法的集 合。GMRES会在前面所有搜索方向上最小化残差,直到重新开始。FGMRE是.GMRES的一个灵活的变种,它能有效地处理更多类的预处理器,通常比GMRES开销2倍多的内存。
3. 常用时间求解器
对于无时间求导的PDE也就是稳态问题,时间求解器有quasi-static method和 Midpoint static Method。一阶时间求导的PDE,如thermal or diffusion equation,有Backward Euler implicit method,backward difference formular methed(BDF)和Generalized midpoint method,其中BDF方法收敛快精度好,推荐使用。二阶时间求导的PDE方法很多,比较代表性的是NEWMark method和HHT alpha method。当然,也可以编写自己定义的求解器,好的求解器衡量标准就是准确度高、速度快。
4. newton求解器
非线性问题都是需要他的,基本方法很简单,就是高等数学第二章的牛顿-莱布尼茨公式:x_new=x+f (x)/f’(x),这种方法属于迭代法,大家都知道迭代法根本思想就是||x_new-x||是否小于某个很小的值,小于就收敛,反之不然。扩展到有限元,f(x)就复杂了一些,成了residual(不要说你不知道什么是residual哦),x成了独立变量。常用的判据有三种:residual增量,x增量,能量增量。其中,能量增量(也就是residual增量乘以x增量)和residual增量是有限元软件的偏好。还有一些变种的newton求解器,原理差不多,修修边幅而已。
这里只是简单提及了常用的求解器,具体的理论和方法在网上和教材上都有,就不赘述了。想要说的是,求解器是数学和计算机背景人非常关心的问题,而放到的力学或者工程应用领域却被极大的忽视了,方便了使用者,但这也许是商用有限元黑盒子带来的弊端。上面有很多求解器都是开源免费的,如果自己编写程序,可以直接从网上下载并使用。
基本上有限元数值方法的最核心内容已经差不多说完了,对应不同的问题人们提出了很多很好的方法,但这些方法都是基于我之前所述的理论基础之上的。接下来谈谈我对有限元将来发展趋势的看法,或许能和打算或正在从事有限元的朋友有些共鸣。
1. 并行计算
随着多核cpu进入个人电脑市场,并行软件已经是大势所趋,那么并行fem软件作为大型科学计算软件,并行趋势极为明朗,各大fem公司分分研制推出并行版,但是并行版的计算准确性与效率还需要经受考验,由于并行机制的设定,每个cpu返回的计算结果不同步,会导致误差。矩阵的blocks分解如何建立行之有效的ghost node/halos,减少communication time,对于各种并行构架的适应性,都是大型有限元并行化要解决的问题,不过当Nvidia公司推出基于Fermi构架的GPU以后,科学计算界为之兴奋,但其是否能催生出有效的并行通用有限元程序还是一个挑战。而且目前个人电脑上并行计算的效率并不能提高很多,有2-3倍就已经很不错了。所以,并行有限元仍需经历考验,但一定是大势所趋。
2. 多尺度模型计算
高性能显微设备的诞生,人们对于微观尺度的物质越发感兴趣,有限元能否描述微观世界呢?纳米级材料,微流场,原子的diffusion,dislocation等等,量子力学由于考虑的电子的作用,使得计算量过于庞大,分子动力学方法虽然已获得了一定成功,但是终归不如有限元来的有效和方便。描述不同尺度下的材料,俨然已成为众多学者们研究的主要方向,在大Z的第6版有限元书中增加了此方面内容,也有不少学者结合向位场,分子动力学和Monte Carlo的方法,来建立多尺度仿真系统,都有一些进展。
后记:我能为有限元做什么
在simwe泡久了,很多朋友发贴子就问,有限元能仿真xx吗?确实,我们都很好奇有限元能做什么,有牛人说给我个棒子我能撬动地球,也有大师说给我一台电脑可以仿真整个世界。但曾几何时我们能为有限元做些什么呢,有的大牛们创造了理论,编写了代码,也有simwe的前辈们建立了这个坛子,编纂了光盘, 作为我自己既不是大师也不是天才少年,只是在有限元里摸爬滚打了几年的老兵,泡泡simwe坛子,回 回帖子,如果本文的某字某句能对您有所帮助,我也就因此感到无比的欣慰了,我想这就是我能为中国有限元事业做得绵薄之事,从而促成了此文。
最后,祝愿所有学习有限元的朋友们,都快乐的生活着。