有限元自动生成及其应用演示

1、1一、通用有限元软件的局限性 通用软件并不通用，只能求解很少的一部分有限元问题。需要有一支训练有素、组织严密、长期稳定的有限元软件队伍，不断维护和改进其通用有限元程序。只能求解今天已有的微分方程，不能求解明天的有限元问题,只能采用今天已出现的算法，不可能采用明天才有的算法。不公开源代码，对用户来说程序是不可改动的，极大地限制了用户的创造性。2二、自己编程的种种困难有限元程序复杂，代码庞大。编程者缺乏软件专业训练，代码可读性差，再用性差，程序难以维护，程序缺乏严密组织与发展。大多数工程师、教授与研究人员，他们主要是研究有限元方法或把有限元方法应用于工程设计与科学研究，而不是致力于有限元软件的开发

2、，因此他们不可能花很多时间编写程序。3三、开发FEPG的动机与目的许多专家希望有自己的有限元程序，但不愿花许多时间编程。FEPG应该不受领域的限制，各种专业的有限元问题都可采用本系统求解。本系统不可能象通用有限元软件那样，通过不断增加单元模块和算法模块来逐步扩大其求解范围。不管什么样的有限元问题，不管采用什么样的算法，都能采用本系统求解。本系统是一种公开源代码的自由软件，任何非本系统开发者都可维护和修改。4四、由有限元公式自动产生有限元程序有限元是求解偏微分方程的数值方法。有限元公式包括： （1）弱形式的微分方程表达式。 （2）形函数表达式。 （3）算法表达式。 a. 形成代数方程组的矩阵表达

3、式。 b. 计算步骤与迭代控制。5五、软件实现思想1、元件化程序设计方法 把有限元程序分解为五个元件程序。6数据转换程序可视化程序后处理程序前处理程序后处理计算程序单元计算程序单元子程序代数方程组求解程序有限元计算部分调用有限元程序结构7有限元程序流程框图线性稳态问题流程前处理程序START元件程序E元件程序SOLV求解器U元件程序DEFER元件程序计算结果图形显示8有限元程序流程框图 线性动态问题流程前处理程序START元件程序BFT元件程序E元件程序SOLV求解器U元件程序DEFER元件程序时间计算是否结束计算结果图形显示否 是9有限元程序流程框图非线性稳态问题流程前处理程序START元件

4、程序E元件程序SOLV求解器U元件程序迭代是否结束计算结果图形显示否 是DEFER元件程序10有限元程序流程框图非线性动态问题流程前处理程序START元件程序BFT元件程序E元件程序SOLV求解器U元件程序DEFER元件程序迭代是否结束计算结果图形显示否 是时间计算是否结束 是否11自动生成系统由近百个元件程序组成。元件程序是以文件名为输入输出参数的可独立运行的程序。由脚本（批命令）文件把元件程序组成一个“程序系统”。大大降低程序的复杂性，大大提高程序的可读性和再用性。实现软件自由化（即自由软件）的一种软件设计方法。12提供一种用以描述有限元问题的语言，以下称为有限元（描述）语言。有限元语言是

5、一种有层次的脚本语言。有限元语言应包括微分方程表达式，形函数表达式，算法表达式等有限元公式的书写形式。有限元语言应与有限元的专业语言尽可能一致，以便广大科技工作者易于理解与接受。有限元语言要适用于各种领域和各种有限元方法与有限元问题。关键问题是如何描述有限元问题，采用什么样的方式准确表达有限元问题及其求解算法。2、人工智能技术13用FEPG生成求解二维线弹性平面应变问题的有限元程序六、程序生成实例14问题的方程描述平衡方程几何方程本构方程15虚功方程虚功方程弱形式将本构关系代入上式即可得以位移为基本未知量的虚功方程16需填写的FEPG公式文件GIO文件VDE文件GCN文件17Strn.gcn文

6、件defia ell &b str astartsin asolvsin astress b18Strn.gio文件dispstrs#elemtype t32dxy19Disp.vde文件disp u vcoor x yfunc ex ey exyshap %1 %2gaus %3mate pe pv fu fv 1.0e10;0.3;0.0;-100.0;$c6 fact = pe/(1.+pv)/(1.-2.*pv)vect u u vvect ed ex eyvect fs fu fvmatr de 2 2(1.-pv) pvpv (1.-pv)20Disp.vde文件funcex=+u

7、/xey=+v/yexy=+u/y+v/xstifdist=+ed_i;ed_j*de_i_j*fact+exy;exy*(0.5-pv)*factload=+u_i*fs_i end21Disp.vde文件funcex=+u/xey=+v/yexy=+u/y+v/xstifdist=+ed_i;ed_j*de_i_j*fact+exy;exy*(0.5-pv)*factload=+u_i*fs_i end22Strs.vde文件disp sa,sb,sc,coef u vcoor x yshap %1 %2gaus %3mass %1mate pe pv fu fv 1.0e10;0.3;0

8、.0;-100.0;$c6 fact = pe/(1.+pv)/(1.-2.*pv)vect u u vvect ed ex eyvect fs fsa fsbvect ss sa sb scvect ff fsa fsb fscmatr de 2 2(1.-pv) pvpv (1.-pv)23Strs.vde文件stif$cv ex=+u/x$cv ey=+v/y$cv fs_i=+de_i_j*ed_j*fact$cv func=+u/y+v/x$c6 fsc=func*(0.5-pv)*factdist=+sa;sa*0.0load=+ss_i*ff_iend24三维正交各向异性材料计算

9、FEPG系统在复合材料中的应用25三维弹性本构关系26其中27其中28需填写的FEPG公式文件GIO文件VDE文件GCN文件29GCN文件defia ell &b str astartsin asolvsin astress b30Gio文件dispstrs#elemtype w43dxyz31Disp.vde文件disp u v wcoor x y zfunc gex gey gez geyz gexz gexyshap %1 %2gaus %3array re(6,6)mate e1 e2 e3 v12 v13 v23 g23 g13 g12 cita fx fy fz 127.6d3;1

10、0.3d3;10.3d3;0.32;0.32;0.32;6.0d3;6.0d3;6.0d3;45;0.0;0.0;5.e3;$c6 pi=3.14159$c6 s=sin(cita*pi/180.0)$c6 c=cos(cita*pi/180.0)$c6 c2=cos(2.0*cita*pi/180.0)$c6 s2=sin(2.0*cita*pi/180.0)fvect le 6fvect ge 6vect gev gex gey gez geyz gexz gexyvect u u v wvect fu fx fy fz32Disp.vde文件matrix te 6 6c*c s*s 0

11、0 0 s2s*s c*c 0 0 0 (-s2/2)0 0 0 1 0 0 0 0 0 c (-s) 00 0 0 s c 0(-s2/2) s2/2 0 0 0 c2matr Re 6 61/e1 (-v12/e1) (-v13/e1) 0 0 0(-v12/e1) 1/e2 (-v23/e2) 0 0 0(-v13/e1) (-v23/e2) 1/e3 0 0 00 0 0 1/g23 0 00 0 0 0 1/g13 00 0 0 0 0 1/g1233Disp.vde文件funcl deform.xyz f le x y z u v wa ge_i=+le_j*te_i_jw gev

12、 gel revmatri n 6 re stifdist=+gev_i;gev_j*re_i_jload=+u_i*fu_iend34平面二维水沙运动问题FEPG系统在水力学中的应用35问题描述 如右图所示，水在1/4环形水渠内流动，进口处x方向单宽流量为4, y方向单宽流量为0,水深为2，出口处水面高度为4。求水渠内的单宽流量和水边高度。 36该问题的偏微分方程水流连续性方程：水流动量方程：37该问题的偏微分方程悬移质不平衡输移方程 悬移质引起的河床变形方程 38算子分裂法求解该问题采用算子分裂法分五步骤求解上述平面二维水沙运动方程 1解水流扩散方程和水位 39需填写的FEPG公式文件GI

13、O文件VDE文件GCN文件40GCN文件defia shola &b sholb c sholcd sholde sholestartsin astartsin bstartsin cstartsin dstartsin ecall trans41if exist stop del stop:1bftsolvsin acopy unod unodaif exist end del end:2solvsin bif not exist end goto 2solvsin csolvsin dsolvsin ecall postif not exist stop goto 142GIO文件fldf

14、lcsddsdczbh#elemtype q42dxy 43FLD.VDE文件disp hu,hv,zcoor x,ycoef hun,hvn,zn,hn,func divvect hu hu hvvect x x yshap %1 %2gaus %3mass %1 volmate epsilon eg en fx fy 0.1d-1 9.8d0 2.0d-2 0.0d0 0.0d0vect fl fx fy44FLD.VDE文件func$c6 vol = 1.0d0$c6 const=1.0div = +hu/x+hv/ystif$cv u2v2=dsqrt(hun*hun+hvn*hvn)

15、dist=+hu_i/x_j;hu_i/x_j*vol*epsilon+div;z*vol+z/x_i;hu_i*vol*eg*hn+hu_i;hu_i*eg*vol*en*en*u2v2/hn*(7.0d0/3.0d0)+z/x_i;z/x_i*vol*constload=+hu_i*vol*fl_iend 45算子分裂法求解该问题2用基于牛顿法的最小二乘法求解水流对流方程 46FLC.VDE文件disp hu hvcoor x y func fhu fhvcoef hun1 hvn1 hun hvn un1 vn1 un vn hn1mate rou 1.0 shap %1 %2gaus

16、%3mass %1 1.0vect hun hun hvn vect x x y vect fhun1 fhun1 fhvn1vect un un vnvect un1 un1 vn1vect hu hu hvvect fhu fhu fhv47func$cv fhun1=+un_j/x_j*hun+hun/x_j*un_j$cv fhun1=+fhun1*dt+hun1$cv fhvn1=+un_j/x_j*hvn+hvn/x_j*un_j$cv fhvn1=+fhvn1*dt+hvn1fhu=+hu+hu_j*un/x_j*dt+hu_j/x_j*un*dt +hu*un_j/x_j*dt

17、+hu/x_j*un_j*dtfhv=+hv+hu_j*vn/x_j*dt+hu_j/x_j*vn*dt +hv*un_j/x_j*dt+hv/x_j*un_j*dtstifdist=+fhu_i;fhu_iload=+fhu_i*fhun1_iend48算子分裂法求解该问题3解挟沙方程扩散项 49SDD.VDE文件disp hscoor x,ycoef hn1,u,vfunc gradhsx gradhsyvect gradhs gradhsx gradhsyvect x x yshap %1 %2gaus %3mass %1 volmate edt rho rous rfai omeg e

18、k eg 0.1d-1 2650 1.0 1.0 1.5d-4 2.2d-1 9.8 50func$c6 vol = 1.0d0gradhsx = +hs/xgradhsy = +hs/y stif$c6 u2v2=dsqrt(u*u+v*v)$c6 ao=rfai*omega$c6 sstar=ek*(u2v2*3.0d0/(eg*hn1*omega)*emdist=+gradhs_i;gradhs_i*vol*edt+hs;hs*rfai*omeg*rho/hn1/rous load=+hs*sstar*rfai*omeg*rho/rous end 51算子分裂法求解该问题4解挟沙方程对流

19、项 52SDC.VDE文件disp hscoor x y func fhs coef hs1 un1 vn1 hn1 hnshap %1 %2gaus %3mate rho 1.0 vect x x y vect un1 un1 vn1 53funcfhs=+hs+hs/x_j*un1_j*dt +hs*un1_j/x_j*dtstifdist=+fhs;fhsload=+fhs*(hs1-rho*(hn-hn1)end 54算子分裂法求解该问题5解河床变形方程 55ZBH.VDE文件disp zbcoor x,ycoef hs,u,v,hshap %1 %2gaus %3mass %1 rh

20、osmate rhos rfai omeg ek eg em 1.4d1 1.0d0 1.5d-4 2.2d-1 9.8d0 7.6d-156stif$c6 u2v2=dsqrt(u*u+v*v)$c6 ao=rfai*omeg$cv sstar=ek*(u2v2*3.0d0/(eg*h*omeg)*emdist=+zb;zb*0.0d0load=+zb*ao*(hs/h-sstar)end 57结果图58线性正弦稳态涡流问题（TEAM Workshop 问题7） FEPG系统在电磁学中的应用59问题描述 TEAM Workshop 问题7是用来检验三维线性正弦涡流场分析方法计算精度的实验模型

21、，其结构如右图所示。该模型由带有方孔的铝板和载流线圈组成，线圈中通有随时间按正弦规律变化的电流，其幅值为2742安匝。铝板的电导率为 求涡流和磁通密度分布。 60求解方程及边界条件61需填写的FEPG公式文件GIO文件VDE文件GCN文件62GCN文件defia ell &b str a &c str a &startc asolvc astress bstress c63GIO文件afimagemagb64方程的afi.vde文件disp Ax,Ay,Az,ucoor x,y,zfunc curlx,curly,curlz,divshap %1 %2gaus %3mass %1 1.0d0m

22、ate sigma epsilon omega fmu fjx fjy fjz 0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;vect A Ax Ay Azvect x x y zvect fj fjx fjy fjzvect curl curlx curly curlz65方程的afi.vde文件disp Ax,Ay,Az,ucoor x,y,zfunc curlx,curly,curlz,divshap %1 %2gaus %3mass %1 1.0d0mate sigma epsilon omega fmu fjx fjy fjz 0.0d0;

23、0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;vect A Ax Ay Azvect x x y zvect fj fjx fjy fjzvect curl curlx curly curlz66afi.vde文件funccurlx=+Az/y-Ay/zcurly=+Ax/z-Az/xcurlz=+Ay/x-Ax/ydiv=+Ax/x+Ay/y+Az/z67afi.vde文件stif$c6 eu1=omega*sigma$c6 eu2=1/fmudist=+curl_i;curl_i*|eu2;0.0|+div;div*|eu2;0.0|+A_i;A_i*|0.

24、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|end68计算涡流的mage.pde文件disp exr exi eyr eyi ezr ezicoef axr axi ayr ayi azr azi ur uicoor x y zshap %1 %2gaus %3mass %1load =fexr fexi feyr feyi fezr fezimate sigma epsilon omega fmu fjx fjy fjz 0.0d0;0.0d0;

25、50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;69mage.pde文件stif$cv fexr = +axi*omega*sigma-ur/x*sigma$cv fexi = -axr*omega*sigma-ui/x*sigma$cv feyr = +ayi*omega*sigma-ur/y*sigma$cv feyi = -ayr*omega*sigma-ui/y*sigma$cv fezr = +azi*omega*sigma-ur/z*sigma$cv fezi = -azr*omega*sigma-ui/z*sigmadist=+exr;exr*0.0end70计算磁通密度的magb.pde文件disp bxr bxi byr byi bzr bzicoef axr axi ayr ayi azr azi ur uicoor x y zshap %1 %2gaus %3mass %1load =fbxr fbxi fbyr fbyi fbzr fbzimate sigma epsilon omega fmu fjx fjy fjz 0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;71magb.pde文件stif$cv fbxr = +a