摘要：多数有限元软件一般只适用于某些特定领域和类型的问题，不能满足其他领域尤其是科研和教学领域工作者的需要，因此开发借助有限元语言等技术自动生成有限元程序的有限元分析软件平台FEPG.介绍用于描述有限元问题的有限元语言：该语言分为两部分，一部分用于描述偏微分方程表达式，另一部分用于描述求解微分方程的有限元算法；用户只需通过有限元语言输入所需的各种表达式和公式即可用FEPG自动产生所需的全部有限元计算源程序；用5个算例说明通过有限元语言编写偏微分方程及其算法的方法.介绍FEPG的生成器技术、组件化技术及其独有的公式库技术；与传统有限元软件采用的程序库技术相比，公式库技术便于用户阅读、修改和维护，能降低软件的复杂度、提高软件的可扩展性.与国内外主流通用有限元软件相比，利用FEPG进行开发还具有程序代码精炼以及投资少等优点，但存在界面不够友好、对用户的起点要求较高等缺点.
关键词：FEPG；有限元；有限元语言；软件开发；软件平台；生成器；组件化；公式库
中图分类号：TB115 文献标志码：B
Finite element analysis software platform FEPG
LIANG Guoping1, TANG Juzhen2
(1. Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China;
2. Beijing FEGEN Software Co., Ltd., Beijing 100190, China)
Abstract： As to the fact that the most finite element softwares are generally applied to certain fields and types of problems, which cannot meet the requirement of workers in the other fields, especially in the fields of research and teaching, finite element analysis software platform FEPG is developed which can automatically generate finite element programs by the techniques such as finite element language. Finite element language is introduced which are used to describe finite element problems. By the language, expressions of partial differential equations and finite element algorithms which are used to solve differential equations can be described; by simply inputting the desired expressions and formulas in finite element language, users can get all the finite element computation source programs by automatic generation function of FEPG; five examples are indicated to introduce the way of preparing partial differential equations and algorithms in finite element language. The main techniques used in FEPG are introduced, including generator, componentization, formula library; Compared with the library of traditional finite element software, formula library is easy to be read, modified and maintained by users and can decrease software complexity and improve software extensibility. Compared with the domestic and foreign major generalpurpose finite element softwares, the development with FEPG can be achieved with less code and investment and short development cycle, but the interfaces of FEPG are not friendly and users are required with a higher starting point.
Key words： FEPG; finite element; finite element language; software development; software platform; generator; componentization; formula library
0引言
随着科学技术和生产力的迅速发展，有限元法的应用越来越广泛，各工矿企业及设计、科研部门已普遍采用有限元法进行生产过程的数值模拟、科研的数值试验和产品的优化设计，有限元法和有限元软件已成为许多高新学科和技术的基本工具和手段，是现代高科技大厦不可缺少的重要支柱之一.
目前国内外许多有限元软件，包括通用有限元程序系统，一般只适用于解决某些特定领域和特定类型的有限元问题，大多根据结构力学的需要设计，对于其他领域的工作者，尤其对于从事科研和教学的工作者，远远不能满足需要.
有限元分析软件平台FEPG的设计思想主要采用组件化的程序设计方法和人工智能技术，根据有限元法的统一数学原理［1］及其内在规律，类似于数学公式的推理方法，根据物理模型和有限元算法，由计算机自动产生有限元程序.
①北京飞箭软件有限公司2003年发布的“FEPG参考手册：基于FEPG的有限元方法”.在FEPG中采用有限元语言［2］求解有限元问题，编写有限元程序的主要工作就是书写有限元问题所需的各种数学表达式，如微分方程表达式、形函数表达式和算法表达式等①；然后，通过FEPG由计算机产生所需的FORTRAN源程序（或其他高级语言的源程序）.因此，有限元语言比目前的高级语言更高一个层次，用它来编写有限元软件能提高效率，并且易于阅读和修改.有限元语言和FEPG从根本上改变有限元软件技术，把广大工程师和科学家从繁琐、重复的有限元编程劳动中彻底解放出来，也十分有利于有限元法在各个领域的推广和应用.
1有限元语言
有限元语言作为一种描述性的脚本模型语言主要用于描述有限元问题，非常接近于有限元书刊论述有限元法和有限元问题时所采用的专业表述语言，简练清晰.
有限元语言分为两部分：一部分用于描述偏微分方程表达式，为此FEPG要求用户按照虚位移原理的形式（弱形式）写一个Partial Differential Equations(PDE)文件，可以采用张量形式，非常符合有限元专业人员的习惯；另一部分用于描述求解微分方程的有限元算法，即描述如何将一个非线性的偏微分方程线性化，如何针对时间依赖问题进行时间离散以及如何控制计算流程和迭代.为此，FEPG要求用户填写两种类型文件：第一种类型文件主要给出线性代数方程组的矩阵和右端项表达式以及迭代步的控制等内容；第二种类型文件给出多物理场采用的算法、耦合方式和有限元计算的流程等.
有限元语言既适用于线性和非线性问题，又适用于稳态和非稳态问题，耦合场的数目可有任意多个，每个场也可有任意多个自由度（即偏微分方程的数目），适用于各种领域和各种工程和学科的有限元问题.
以下通过5个算例说明如何通过有限元语言书写偏微分方程及其算法.
算例1考虑稳态二维热传导问题
泊松方程为-•(ku)=qinΩ边界条件为u(x,y)=g,onΓ1
-kun=a(u－b)， onΓ2其有限元弱形式为
∫Ωkuxδux+uyδuydΩ=∫ΩqδudΩ+∫ΩkunδudΩ
将上面方程的体积分项写成PDE文件如下：
disp u,
coor x,y,
shap q 4
gaus q
mate ek eq 1.0d0;0.0;
stif
dist =+［u/x;u/x］*ek+［u/y;u/y］*ek
load =+［u］*eq
end
算法描述GCN文件如下：
defi
a ell
Startsin a
Solvsin a
① 北京飞箭软件有限公司2003年发布的“FEPG参考手册：基于FEPG的有限元方法”.算例2线弹性小变形的PDE文件描述
平衡方程为σxxx+σxyy+fx=0
σxyx+σyyy+fy=0几何方程为εxx=ux,εyy=vy,εxy=uy+vx本构方程为
σxx
σyy
σxy=E(1+ν)(1－2ν)•
1－νν0
ν1－ν0
000.5－νεxx
εyy
εxy
第一类边界条件为u=u0，v=v0第二类边界条件为Tx=f1，Ty=f2第三类边界条件为Tx=f1(u,v,w)，Ty=f2(u,v,w)运用迦辽金有限元法求位移，由平衡方程得
∫Vσxxx+σxyy+fxδu+
σxyx+σyyy+fyδvdV=0
将其转化为弱形式，得
∫ V(σxxδεxx+σyyδεyy+σxyδεxy)dV=
∫ V(fxδu+fyδv)dV+∫ Γ(Txδu+Tyδv)dΓ
将本构关系代入上式即可得以位移为基本未知量的虚功方程.将上述方程的体积分项写成PDE文件①如下：
disp u v
coor x y
func ex ey exy
shap %1 %2
gaus %3
mate pe pv fu fv
S|c6 fact = pe/(1.+pv)/(1.-2.*pv)
vect u u v
vect ed ex ey
vect fs fu fv
matr de 2 2
(1.-pv) pv
pv (1.-pv)
func
ex=+［u/x］
ey=+［v/y］
exy=+［u/y］+［v/x］
stif
dist=+［ed_i;ed_j］*de_i_j*fact +［exy;exy］*
(0.5-pv)*fact
load=+［u_i］*fs_i
end
算例3对流项采用有限体积法求解流体力学NS方程为)u+p=f
•u=0对流项采用迎风有限体积法，其他项采用有限元法，其PDE文件如下：
defi
disp u v p
coor x y
coef un vn
func div
shap q 9
shap q 4 p 1 2 3 4
gaus q
mate emu rou fx fy ＼1.0d-3 1.0d0 0 0
gvs q4
func
div=+［u/x］+［v/y］
stif
dist=+［u/x;u/x］*emu+［u/y;u/y］*emu+
［v/x;v/x］*emu+［v/y;v/y］*emu+［div;p］-［p;div］+［p/x;p/x］*1.0d-10+
［p/y;p/y］*1.0d-10
load=+［u］*fx*rou+［v］*fy*rou
End
算例4土坝的三维无压渗流-•(Kφ)=Q式中：φ为水头；K为渗透因数，K=1.0；Q为源汇，Q=0.
由于自由渗流面未知，这是一个非线性问题.采用死活节点法求解方程-•(Kφ)+mφ=Q
φ≥y，K=0,m=ε
φ defi disp u coor x y z coef un shap w 4 gaus %3 mate pk eq 1.0 0.0 stif S|c6 pk1=pk S|c6 em=0.0d0 S|c6 if(un.lt.y+1.0d-11) then S|c6 pk1=0.0d0 S|c6 em=1.0d-8 S|c6 endif dist=+［u/x;u/x］*pk1+［u/y;u/y］*pk1+ ［u/z;u/z］*pk1+［u;u］*em load=-［u］*eq end 计算结果见图1. 图 1土坝三维无压渗流场分析 Fig.1Analysis on 3D unconfined seepage flow field of earth dam 算例5三维电磁场A-φ方法的PDE文件 ×1μ×A-1μ•A+ σjωA+σφ=0涡流区 •-σjωA-σφ=0 ×1μ×A-1μ•A=Js非涡流区 边界条件为A=0，φ=0无限远处 n•(-σjωA-σφ)=0 导体和非导体的交界面处PDE文件如下： defi disp Ax Ay Az u coor x y z func curlx curly curlz div shap %1 %2 gaus %3 mass %1 1.0d0 mate sigma epsilon omega fmu fjx fjy fjz vect A Ax Ay Az vect x x y z vect fj fjx fjy fjz vect curl curlx curly curlz func curlx=+［Az/y］-［Ay/z］ curly=+［Ax/z］-［Az/x］ curlz=+［Ay/x］-［Ax/y］ div=+［Ax/x］+［Ay/y］+［Az/z］ stif S|c6 eu1=omega*sigma S|c6 eu2=1/fmu dist=+［curl_i;curl_i］*|eu2;0.0|+［div;div］* |eu2;0.0|+［A_i;A_i］*|0.0;eu1|+ ［u/x_i;A_i］*|sigma;0.0|+［A_i;u/x_i］* |0.0;eu1|+［u/x_i;u/x_i］*|sigma;0.0| load=+［A_i］*fj_i*|1.0;0.0| end 2三项软件技术 FEPG采用生成器技术（见图2）、组件化技术［3］和公式库技术等三项软件技术. 图 2生成器技术 Fig.2Generator technology 采用有限元语言编程就是书写微分方程表达式和算法，很容易阅读、理解、修改和再用，并且需要书写的代码量与由它产生的FORTRAN语言有限元程序相比减少90%以上，但灵活性却相同.由此可见，FEPG可为用户节省大量的编程时间，让用户将主要精力集中在最有创造性的工作、物理模型建立（即偏微分方程表达式）和算法研究上. 组件化是当今许多软件公司为提高软件的再用性而普遍采用的一项软件技术.这项技术的关键是组件之间接口的统一规定，FEPG采用文件名作为组件之间的接口.以文件名为接口书写简单，内容丰富、灵活. 众所周知，许多软件系统（包括通用有限元软件和平台软件）为使功能更加完善、求解的问题更加广泛，都开发出由各种子程序组成的庞大的程序库，但FEPG不同，它主要依靠的不是程序库而是公式库，由各种各样的有限元公式组成的公式库包括形函数库、微分算子库、PDE模型库、有限元算法库以及各学科领域的专业公式库，包括结构力学的杆、梁、板和壳等，固体力学的线弹性、大变形、弹塑性和黏弹性等，流体力学的NS方程、Dacy流和非牛顿流等，电磁场的Aφ方法和棱单元法等. 公式库技术是FEPG独有的技术，与程序库技术相比，它不会因为库功能的扩大而增加软件的复杂性，而且更容易阅读、修改和维护，其可重用性显而易见，并且通过公式库技术，用户更容易参与软件的发展. 为把分别适用于固体力学的有限元法、流体力学的有限体积法和电磁场计算的棱单元等有机地融为一体，并集成到统一的数值仿真程序中，FEPG对这些方法的软件实现进行独创性的改造. 3与通用软件比较 目前，国内外已出现许多功能强大的通用有限元软件，如ADINA，ANSYS，Abaqus，MSC Marc和COMSOL等软件.对于一般的科研人员，这些通用程序仅能解决某一类型或某几种类型的问题，而FEPG可以生成解决任何有限元问题的程序，编程人员可以拥有生成程序的自主版权. 相对于通用有限元软件，FEPG亦存在不足之处，如界面不够友好、对用户的起点要求较高以及难学难用等.因而，有待进一步改善FEPG的用户界面、完善公式库等，使其能兼具通用有限元软件易学易用的优点.FEPG与通用有限元软件的比较见表1. 表 1FEPG与通用有限元软件的比较 Tab.1Comparison of FEPG with general finite element software通用有限元软件FEPG程序代码量庞大（几十万行）程序代码量精炼（几万行）开发投资大开发投资少开发周期长开发周期短源代码不可见，难修改源代码可见，易修改难维护易维护难集成易集成面向“傻瓜”面向科研“自主创新”解决的问题模型单一解决问题模型复杂多变，灵活性强易学易用难学难用 4结论 FEPG采用有限元语言编程，代码量比生成器产生的FORTRAN语言的代码量减少90%以上.采用张量形式表示弱形式的PDE符合有限元专业人员的习惯,易阅读、理解、修改和维护.有限元语言是一种超平台语言，它既不依赖操作系统，也不依赖计算机；既适用于单CPU的串行机，也适用于多CPU的并行机；适用于各学科领域及各种有限元问题；为科学家和工程师提供求解偏微分方程的有力工具. 组件化技术与公式库技术的应用，使得人人都可以在有限元分析软件平台上增加、修改和更换组件程序和公式库. FEPG是第一个能通过互联网使用的有限元软件平台.通过互联网使用，服务器只需生成FORTRAN程序，速度快、代码少、网络传输量少.生成的代码可在任何一台串行或并行的计算机上编译运行. 参考文献： ［1］ZIENKIEWICZ Q C, TAYLOR R L, ZHU J Z. The finite element method: its basis and fundamentals［M］. 6th ed. Oxford: ButterworthHeinemann, 2008: 733. ［2］梁国平. 有限元语言［M］. 北京: 科学出版社, 2009: 180. ［3］梁国平, 傅子智. 程序元件化方法及元件程序组装语言［J］. 计算结构力学及其应用, 1985, 2(2): 121125. LIANG Guoping, FU Zizhi. The component programming method and its language［J］. Comput Struct Mech & Applications, 1985, 2(2): 121125. (编辑于杰)