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

《有限元的自动生成及其应用演示》由会员分享，可在线阅读，更多相关《有限元的自动生成及其应用演示（89页珍藏版）》请在装配图网上搜索。

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

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

3、1）弱形式的微分方程表达式。 （2）形函数表达式。 （3）算法表达式。 a. 形成代数方程组的矩阵表达式。 b. 计算步骤与迭代控制。五、软件实现思想五、软件实现思想1、元件化程序设计方法、元件化程序设计方法 把有限元程序分解为五个元件程序。数据转换程序可视化程序后 处 理程序前处理程序后处理计算程序单元计算程序单元子程序代数方程组求解程序有限元计算部分调用有限元程序结构有限元程序结构有限元程序流程框图有限元程序流程框图线性稳态问题流程前处理程序START元件程序E元件程序SOLV求解器U元件程序DEFER元件程序计算结果图形显示有限元程序流程框图有限元程序流程框图 线性动态问题流程前处理程序

4、START元件程序BFT元件程序E元件程序SOLV求解器U元件程序DEFER元件程序时间计算是否结束计算结果图形显示否 是有限元程序流程框图有限元程序流程框图非线性稳态问题流程前处理程序START元件程序E元件程序SOLV求解器U元件程序迭代是否结束计算结果图形显示否 是DEFER元件程序有限元程序流程框图有限元程序流程框图非线性动态问题流程前处理程序START元件程序BFT元件程序E元件程序SOLV求解器U元件程序DEFER元件程序迭代是否结束计算结果图形显示否 是时间计算是否结束 是否自动生成系统由近百个元件程序组成。元件程序是以文件名为输入输出参数的可独立运行的程序。由脚本（批命令）文件

5、把元件程序组成一个“程序系统”。大大降低程序的复杂性，大大提高程序的可读性和再用性。实现软件自由化（即自由软件）的一种软件设计方法。提供一种用以描述有限元问题的语言，以下称为有限元（描述）语言。有限元语言是一种有层次的脚本语言。有限元语言应包括微分方程表达式，形函数表达式，算法表达式等有限元公式的书写形式。有限元语言应与有限元的专业语言尽可能一致，以便广大科技工作者易于理解与接受。有限元语言要适用于各种领域和各种有限元方法与有限元问题。关键问题是如何描述有限元问题，采用什么样的方式准确表达有限元问题及其求解算法。2、人工智能技术、人工智能技术用用FEPG生成求解二维线弹性生成求解二维线弹性平面

6、应变问题的有限元程序平面应变问题的有限元程序六、程序生成实例六、程序生成实例问题的方程描述平衡方程00yyyxyxxyxxfyxfyxxvyuyvxuxyyyxx,xyyyxxxyyyxxE5 . 0000101)21)(1 (几何方程本构方程虚功方程虚功方程虚功方程弱形式虚功方程弱形式 )()()(dvTuTdVvfufdVyxVyxVxyxyyyyyxxxx0)()( VyyyxyxxyxxdVvfyxufyx将本构关系代入上式即可得以位移为基本未知量的虚功方程需填写的需填写的FEPG公式文件公式文件GIO文件文件VDE文件文件GCN文件文件Strn.gcn文件文件defidefia el

7、l &a ell &b str ab str astartsin astartsin asolvsin asolvsin astress bstress bStrn.gio文件文件dispdispstrsstrs#elemtype t3#elemtype t32dxy2dxyDisp.vde文件文件disp u vdisp u vcoor x ycoor x yfunc ex ey exyfunc ex ey exyshap %1 %2shap %1 %2gaus %3gaus %3mate pe pv fu fv 1.0e10;0.3;0.0;-100.0;mate pe pv fu fv

8、1.0e10;0.3;0.0;-100.0;$c6 fact = pe/(1.+pv)/(1.-2.$c6 fact = pe/(1.+pv)/(1.-2.* *pv)pv)vect u u vvect u u vvect ed ex eyvect ed ex eyvect fs fu fvvect fs fu fvmatr de 2 2matr de 2 2(1.-pv) pv(1.-pv) pvpv (1.-pv)pv (1.-pv)Disp.vde文件文件funcfuncex=+u/xex=+u/xey=+v/yey=+v/yexy=+u/y+v/xexy=+u/y+v/xstifsti

9、fdist=+ed_i;ed_jdist=+ed_i;ed_j* *de_i_jde_i_j* *fact+exy;exyfact+exy;exy* *(0.5-pv)(0.5-pv)* *factfactload=+u_iload=+u_i* *fs_i fs_i endendDisp.vde文件文件funcfuncex=+u/xex=+u/xey=+v/yey=+v/yexy=+u/y+v/xexy=+u/y+v/xstifstifdist=+ed_i;ed_jdist=+ed_i;ed_j* *de_i_jde_i_j* *fact+exy;exyfact+exy;exy* *(0.5-

10、pv)(0.5-pv)* *factfactload=+u_iload=+u_i* *fs_i fs_i endendStrs.vde文件文件disp sa,sb,sc,disp sa,sb,sc,coef u vcoef u vcoor x ycoor x yshap %1 %2shap %1 %2gaus %3gaus %3mass %1mass %1mate pe pv fu fv 1.0e10;0.3;0.0;-100.0;mate pe pv fu fv 1.0e10;0.3;0.0;-100.0;$c6 fact = pe/(1.+pv)/(1.-2.$c6 fact = pe/(

11、1.+pv)/(1.-2.* *pv)pv)vect u u vvect u u vvect ed ex eyvect ed ex eyvect fs fsa fsbvect fs fsa fsbvect ss sa sb scvect ss sa sb scvect ff fsa fsb fscvect ff fsa fsb fscmatr de 2 2matr de 2 2(1.-pv) pv(1.-pv) pvpv (1.-pv)pv (1.-pv)Strs.vde文件文件stifstif$cv ex=+u/x$cv ex=+u/x$cv ey=+v/y$cv ey=+v/y$cv fs

12、_i=+de_i_j$cv fs_i=+de_i_j* *ed_jed_j* *factfact$cv func=+u/y+v/x$cv func=+u/y+v/x$c6 fsc=func$c6 fsc=func* *(0.5-pv)(0.5-pv)* *factfactdist=+sa;sadist=+sa;sa* *0.00.0load=+ss_iload=+ss_i* *ff_iff_iendend三维正交各向异性材料计算三维正交各向异性材料计算FEPG系统在复合材料中的应用三维弹性本构关系三维弹性本构关系xyxzyzzyxxyxyxzyzzyxDTxyTDTD其中其中222222000

13、0000000000010020002000nmmnmnmnnmmnmnmnnmTsin,cosnm其中其中1232323321221111/1000000/1000000/1000000/1/23/13000/23/1/12000/13/12/1GGGEEvEvEvEEvEvEvED需填写的需填写的FEPG公式文件公式文件GIO文件文件VDE文件文件GCN文件文件GCN文件文件defidefia ell &a ell &b str ab str astartsin astartsin asolvsin asolvsin astress bstress bGio文件文件dispdispstrs

14、strs#elemtype w4#elemtype w43dxyz3dxyzDisp.vde文件文件disp u v wdisp u v wcoor x y zcoor x y zfunc gex gey gez geyz gexz gexyfunc gex gey gez geyz gexz gexyshap %1 %2shap %1 %2gaus %3gaus %3array re(6,6)array re(6,6)mate e1 e2 e3 v12 v13 v23 g23 g13 g12 cita fx fy fz mate e1 e2 e3 v12 v13 v23 g23 g13 g1

15、2 cita fx fy fz 127.6d3;10.3d3;10.3d3;0.32;0.32;0.32;6.0d3;6.0d3;6.0d3;45127.6d3;10.3d3;10.3d3;0.32;0.32;0.32;6.0d3;6.0d3;6.0d3;45;0.0;0.0;5.e3;0.0;0.0;5.e3;$c6 pi=3.14159$c6 pi=3.14159$c6 s=sin(cita$c6 s=sin(cita* *pi/180.0)pi/180.0)$c6 c=cos(cita$c6 c=cos(cita* *pi/180.0)pi/180.0)$c6 c2=cos(2.0$c6

16、 c2=cos(2.0* *citacita* *pi/180.0)pi/180.0)$c6 s2=sin(2.0$c6 s2=sin(2.0* *citacita* *pi/180.0)pi/180.0)fvect le 6fvect le 6fvect ge 6fvect ge 6vect gev gex gey gez geyz gexz gexyvect gev gex gey gez geyz gexz gexyvect u u v wvect u u v wvect fu fx fy fzvect fu fx fy fzDisp.vde文件文件matrix te 6 6matrix

17、 te 6 6c c* *c sc s* *s 0 0 0 s2s 0 0 0 s2s s* *s cs c* *c 0 0 0 (-s2/2)c 0 0 0 (-s2/2)0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 c (-s) 00 0 0 c (-s) 00 0 0 s c 00 0 0 s c 0(-s2/2) s2/2 0 0 0 c2(-s2/2) s2/2 0 0 0 c2matr Re 6 6matr Re 6 61/e1 (-v12/e1) (-v13/e1) 0 0 01/e1 (-v12/e1) (-v13/e1) 0 0 0(-v12/e1) 1/e2

18、(-v23/e2) 0 0 0(-v12/e1) 1/e2 (-v23/e2) 0 0 0(-v13/e1) (-v23/e2) 1/e3 0 0 0(-v13/e1) (-v23/e2) 1/e3 0 0 00 0 0 1/g23 0 00 0 0 1/g23 0 00 0 0 0 1/g13 00 0 0 0 1/g13 00 0 0 0 0 1/g120 0 0 0 0 1/g12Disp.vde文件文件funcfuncl deform.xyz f le x y z u v wl deform.xyz f le x y z u v wa ge_i=+le_ja ge_i=+le_j* *

19、te_i_jte_i_jw gev gew gev gel revmatri n 6 rel revmatri n 6 re stifstifdist=+gev_i;gev_jdist=+gev_i;gev_j* *re_i_jre_i_jload=+u_iload=+u_i* *fu_ifu_iendend平面二维水沙运动问题平面二维水沙运动问题FEPG系统在水力学中的应用问题描述问题描述 如右图所示，水在1/4环形水渠内流动，进口处x方向单宽流量为4, y方向单宽流量为0,水深为2，出口处水面高度为4。求水渠内的单宽流量和水边高度。 该问题的偏微分方程该问题的偏微分方程0)()(yHVxH

20、UtZxfHVUUgnxZgHyHUyxHUxyHUVxHUUtHU3/ 1222)()()()()(yfHVUVgnyZgHyHVyxHVxyHVVxHVUtHV3/1222)()()()()(水流连续性方程：水流连续性方程：水流动量方程：水流动量方程：该问题的偏微分方程该问题的偏微分方程悬移质不平衡输移方程悬移质不平衡输移方程 悬移质引起的河床变形方程悬移质引起的河床变形方程 yHSDyxHSDxtZyHSVxHSUtHStt)()()()()(*SStZbsmgHVUkS322*算子分裂法求解该问题算子分裂法求解该问题采用算子分裂法分五步骤求解上述平面二维水沙运动方程 1解水流扩散方程和

21、水位解水流扩散方程和水位 iijijiiifxgHxwxtwxwt0需填写的需填写的FEPG公式文件公式文件GIO文件文件VDE文件文件GCN文件文件GCN文件文件defidefia shola &a shola &b sholb b sholb c sholcc sholcd sholdd sholde sholee shole startsin astartsin astartsin bstartsin bstartsin cstartsin cstartsin dstartsin dstartsin estartsin ecall transcall transif exist stop

22、 del stopif exist stop del stop:1:1bftbftsolvsin asolvsin acopy unod unodacopy unod unodaif exist end del endif exist end del end:2:2solvsin bsolvsin bif not exist end goto 2if not exist end goto 2solvsin csolvsin csolvsin dsolvsin dsolvsin esolvsin ecall postcall postif not exist stop goto 1if not

23、exist stop goto 1GIO文件文件fldfldflcflcsddsddsdcsdczbhzbh#elemtype q4#elemtype q42dxy 2dxy FLD.VDE文件文件disp hu,hv,zdisp hu,hv,zcoor x,ycoor x,ycoef hun,hvn,zn,hn,coef hun,hvn,zn,hn,func divfunc divvect hu hu hvvect hu hu hvvect x x yvect x x yshap %1 %2shap %1 %2gaus %3gaus %3mass %1 volmass %1 volmate

24、epsilon eg en fx fy 0.1d-1 9.8d0 2.0d-2 0.0d0 0.0d0mate epsilon eg en fx fy 0.1d-1 9.8d0 2.0d-2 0.0d0 0.0d0vect fl fx fyvect fl fx fyFLD.VDE文件文件funcfunc$c6 vol = 1.0d0$c6 vol = 1.0d0$c6 const=1.0$c6 const=1.0div = +hu/x+hv/ydiv = +hu/x+hv/y stifstif$cv u2v2=dsqrt(hun$cv u2v2=dsqrt(hun* *hun+hvnhun+h

25、vn* *hvn)hvn)dist=+hu_i/x_j;hu_i/x_jdist=+hu_i/x_j;hu_i/x_j* *volvol* *epsilonepsilon+div;z+div;z* *vol+z/x_i;hu_ivol+z/x_i;hu_i* *volvol* *egeg* *hnhn+hu_i;hu_i+hu_i;hu_i* *egeg* *volvol* *enen* *enen* *u2v2/hnu2v2/hn* * *(7.0d0/3.0d0)(7.0d0/3.0d0)+z/x_i;z/x_i+z/x_i;z/x_i* *volvol* *constconst load

26、=+hu_iload=+hu_i* *volvol* *fl_ifl_i end end 算子分裂法求解该问题算子分裂法求解该问题2用基于牛顿法的最小二乘法求解水流对流方程用基于牛顿法的最小二乘法求解水流对流方程 jijijjnijjjijjjjijjiixwuwxutwxwuwxuxwuwxutwFLC.VDE文件文件disp hu hvdisp hu hvcoor x y coor x y func fhu fhvfunc fhu fhvcoef hun1 hvn1 hun hvn un1 vn1 un vn hn1coef hun1 hvn1 hun hvn un1 vn1 un vn

27、hn1mate rou 1.0 mate rou 1.0 shap %1 %2shap %1 %2gaus %3gaus %3mass %1 1.0mass %1 1.0vect hun hun hvn vect hun hun hvn vect x x y vect x x y vect fhun1 fhun1 fhvn1vect fhun1 fhun1 fhvn1vect un un vnvect un un vnvect un1 un1 vn1vect un1 un1 vn1vect hu hu hvvect hu hu hvvect fhu fhu fhvvect fhu fhu fh

28、vfuncfunc$cv fhun1=+un_j/x_j$cv fhun1=+un_j/x_j* *hun+hun/x_jhun+hun/x_j* *un_jun_j$cv fhun1=+fhun1$cv fhun1=+fhun1* *dt+hun1dt+hun1$cv fhvn1=+un_j/x_j$cv fhvn1=+un_j/x_j* *hvn+hvn/x_jhvn+hvn/x_j* *un_jun_j$cv fhvn1=+fhvn1$cv fhvn1=+fhvn1* *dt+hvn1dt+hvn1fhu=+hu+hu_jfhu=+hu+hu_j* *un/x_jun/x_j* *dt+

29、hu_j/x_jdt+hu_j/x_j* *unun* *dtdt +hu +hu* *un_j/x_jun_j/x_j* *dt+hu/x_jdt+hu/x_j* *un_jun_j* *dtdt fhv=+hv+hu_jfhv=+hv+hu_j* *vn/x_jvn/x_j* *dt+hu_j/x_jdt+hu_j/x_j* *vnvn* *dtdt +hv +hv* *un_j/x_jun_j/x_j* *dt+hv/x_jdt+hv/x_j* *un_jun_j* *dtdt stifstifdist=+fhu_i;fhu_idist=+fhu_i;fhu_i load=+fhu_il

30、oad=+fhu_i* *fhun1_ifhun1_i endend算子分裂法求解该问题算子分裂法求解该问题3解挟沙方程扩散项解挟沙方程扩散项 itiSxHSDxSSttHS*SDD.VDE文件文件disp hsdisp hscoor x,ycoor x,ycoef hn1,u,vcoef hn1,u,vfunc gradhsx gradhsyfunc gradhsx gradhsyvect gradhs gradhsx gradhsyvect gradhs gradhsx gradhsyvect x x yvect x x yshap %1 %2shap %1 %2gaus %3gaus %

31、3mass %1 volmass %1 volmate edt rho rous rfai omeg ek eg mate edt rho rous rfai omeg ek eg 0.1d-1 2650 1.0 1.0 1.5d-4 2.2d-1 9.8 0.1d-1 2650 1.0 1.0 1.5d-4 2.2d-1 9.8 funcfunc$c6 vol = 1.0d0$c6 vol = 1.0d0gradhsx = +hs/xgradhsx = +hs/x gradhsy = +hs/ygradhsy = +hs/y stifstif$c6 u2v2=dsqrt(u$c6 u2v2=

32、dsqrt(u* *u+vu+v* *v)v)$c6 ao=rfai$c6 ao=rfai* *omegaomega$c6 sstar=ek$c6 sstar=ek* *(u2v2(u2v2* * *3.0d0/(eg3.0d0/(eg* *hn1hn1* *omega)omega)* * *ememdist=+gradhs_i;gradhs_idist=+gradhs_i;gradhs_i* *volvol* *edt+hs;hsedt+hs;hs* *rfairfai* *omegomeg* *rho/hn1/rho/hn1/rousrous load=+hsload=+hs* *ssta

33、rsstar* *rfairfai* *omegomeg* *rho/rousrho/rous end end 算子分裂法求解该问题算子分裂法求解该问题4解挟沙方程对流项解挟沙方程对流项 0iixHSutHtHSSDC.VDE文件文件disp hsdisp hscoor x y coor x y func fhs func fhs coef hs1 un1 vn1 hn1 hncoef hs1 un1 vn1 hn1 hnshap %1 %2shap %1 %2gaus %3gaus %3mate rho 1.0 mate rho 1.0 vect x x y vect x x y vect

34、un1 un1 vn1 vect un1 un1 vn1 funcfuncfhs=+hs+hs/x_jfhs=+hs+hs/x_j* *un1_jun1_j* *dtdt +hs +hs* *un1_j/x_jun1_j/x_j* *dtdt stifstifdist=+fhs;fhsdist=+fhs;fhs load=+fhsload=+fhs* *(hs1-rho(hs1-rho* *(hn-hn1)(hn-hn1) end end 算子分裂法求解该问题算子分裂法求解该问题5解河床变形方程解河床变形方程 )(*SStZbsZBH.VDE文件文件disp zbdisp zbcoor x,y

35、coor x,ycoef hs,u,v,hcoef hs,u,v,hshap %1 %2shap %1 %2gaus %3gaus %3mass %1 rhosmass %1 rhosmate rhos rfai omeg ek eg em mate rhos rfai omeg ek eg em 1.4d1 1.0d0 1.5d-4 2.2d-1 9.8d0 7.6d-11.4d1 1.0d0 1.5d-4 2.2d-1 9.8d0 7.6d-1 stifstif$c6 u2v2=dsqrt(u$c6 u2v2=dsqrt(u* *u+vu+v* *v)v)$c6 ao=rfai$c6 ao

36、=rfai* *omegomeg$cv sstar=ek$cv sstar=ek* *(u2v2(u2v2* * *3.0d0/(eg3.0d0/(eg* *h h* *omeg)omeg)* * *ememdist=+zb;zbdist=+zb;zb* *0.0d00.0d0 load=+zbload=+zb* *aoao* *(hs/h-sstar)(hs/h-sstar) end end 结果图结果图线性正弦稳态涡流问题线性正弦稳态涡流问题（TEAM Workshop TEAM Workshop 问题问题7 7） FEPG系统在电磁学中的应用问题描述问题描述 TEAM WorkshopT

37、EAM Workshop 问题7是用来检验三维线性正弦涡流场分析方法计算精度的实验模型，其结构如右图所示。该模型由带有方孔的铝板和载流线圈组成，线圈中通有随时间按正弦规律变化的电流，其幅值为2742安匝。铝板的电导率为 求涡流和磁通密度分布。 713.526 10 ()m求解方程求解方程及边界条件及边界条件110011sjjJ 方程：在涡流区在非涡流区000nj 边界条件：，在无限远处在导体和非导体的交界面处需填写的需填写的FEPG公式文件公式文件GIO文件文件VDE文件文件GCN文件文件GCN文件文件defidefia ell &a ell &b str a &b str a &c str

38、a &c str a &startc astartc asolvc asolvc astress bstress bstress cstress cGIO文件文件afiafimagemagemagbmagb方程的方程的afi.vde文件文件disp Ax,Ay,Az,udisp Ax,Ay,Az,ucoor x,y,zcoor x,y,zfunc curlx,curly,curlz,divfunc curlx,curly,curlz,divshap %1 %2shap %1 %2gaus %3gaus %3mass %1 1.0d0mass %1 1.0d0mate sigma epsilon

39、 omega fmu fjx fjy fjz mate sigma epsilon omega fmu fjx fjy fjz 0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;vect A Ax Ay Azvect A Ax Ay Azvect x x y zvect x x y zvect fj fjx fjy fjzvect fj fjx fjy fjzvect curl curlx curly curlzvect curl curlx curly curlz

40、方程的方程的afi.vdeafi.vde文件文件disp Ax,Ay,Az,udisp Ax,Ay,Az,ucoor x,y,zcoor x,y,zfunc curlx,curly,curlz,divfunc curlx,curly,curlz,divshap %1 %2shap %1 %2gaus %3gaus %3mass %1 1.0d0mass %1 1.0d0mate sigma epsilon omega fmu fjx fjy fjz mate sigma epsilon omega fmu fjx fjy fjz 0.0d0;0.0d0;50.d0;1.256d-6;0.0d0

41、;0.0d0;0.0d0;0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;vect A Ax Ay Azvect A Ax Ay Azvect x x y zvect x x y zvect fj fjx fjy fjzvect fj fjx fjy fjzvect curl curlx curly curlzvect curl curlx curly curlzafi.vde文件文件funcfunccurlx=+Az/y-Ay/zcurlx=+Az/y-Ay/zcurly=+Ax/z-Az/xcurly=+Ax/z-Az/xcurlz=+Ay/x-A

42、x/ycurlz=+Ay/x-Ax/ydiv=+Ax/x+Ay/y+Az/zdiv=+Ax/x+Ay/y+Az/zafi.vde文件文件stifstif$c6 eu1=omega$c6 eu1=omega* *sigmasigma$c6 eu2=1/fmu$c6 eu2=1/fmudist=+curl_i;curl_idist=+curl_i;curl_i* *|eu2;0.0|eu2;0.0|+div;div+div;div* *|eu2;0.0|eu2;0.0|+A_i;A_i+A_i;A_i* *|0.0;eu1|0.0;eu1|+u/x_i;A_i+u/x_i;A_i* *|sigma

43、;0.0|sigma;0.0|+A_i;u/x_i+A_i;u/x_i* *|0.0;eu1|0.0;eu1|+u/x_i;u/x_i+u/x_i;u/x_i* *|sigma;0.0|sigma;0.0|load=+A_iload=+A_i* *fj_ifj_i* *|1.0;0.0|1.0;0.0|endend计算涡流的mage.pde文件disp exr exi eyr eyi ezr ezidisp exr exi eyr eyi ezr ezicoef axr axi ayr ayi azr azi ur uicoef axr axi ayr ayi azr azi ur uicoo

44、r x y zcoor x y zshap %1 %2shap %1 %2gaus %3gaus %3mass %1mass %1load =fexr fexi feyr feyi fezr feziload =fexr fexi feyr feyi fezr fezimate sigma epsilon omega fmu fjx fjy fjz mate sigma epsilon omega fmu fjx fjy fjz 0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0

45、;0.0d0;mage.pde文件文件stifstif$cv fexr = +axi$cv fexr = +axi* *omegaomega* *sigma-ur/xsigma-ur/x* *sigmasigma$cv fexi = -axr$cv fexi = -axr* *omegaomega* *sigma-ui/xsigma-ui/x* *sigmasigma$cv feyr = +ayi$cv feyr = +ayi* *omegaomega* *sigma-ur/ysigma-ur/y* *sigmasigma$cv feyi = -ayr$cv feyi = -ayr* *ome

46、gaomega* *sigma-ui/ysigma-ui/y* *sigmasigma$cv fezr = +azi$cv fezr = +azi* *omegaomega* *sigma-ur/zsigma-ur/z* *sigmasigma$cv fezi = -azr$cv fezi = -azr* *omegaomega* *sigma-ui/zsigma-ui/z* *sigmasigmadist=+exr;exrdist=+exr;exr* *0.00.0endend计算磁通密度的计算磁通密度的magb.pde文件文件disp bxr bxi byr byi bzr bzidisp

47、 bxr bxi byr byi bzr bzicoef axr axi ayr ayi azr azi ur uicoef axr axi ayr ayi azr azi ur uicoor x y zcoor x y zshap %1 %2shap %1 %2gaus %3gaus %3mass %1mass %1load =fbxr fbxi fbyr fbyi fbzr fbziload =fbxr fbxi fbyr fbyi fbzr fbzimate sigma epsilon omega fmu fjx fjy fjz mate sigma epsilon omega fmu

48、fjx fjy fjz 0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;magb.pde文件文件stifstif$cv fbxr = +azr/y-ayr/z$cv fbxr = +azr/y-ayr/z$cv fbxi = +azi/y-ayi/z$cv fbxi = +azi/y-ayi/z$cv fbyr = +axr/z-azr/x$cv fbyr = +axr/z-azr/x$cv fbyi = +axi/z-azi/x$cv fbyi = +axi/z

49、-azi/x$cv fbzr = +ayr/x-axr/y$cv fbzr = +ayr/x-axr/y$cv fbzi = +ayi/x-axi/y$cv fbzi = +ayi/x-axi/ydist=+bxr;bxrdist=+bxr;bxr* *0.00.0endend 结果图结果图导体表面导体表面z=0.019mz=0.019m处涡流的矢量图处涡流的矢量图 结果图结果图BzBz的计算值与测量值的比较图（的计算值与测量值的比较图（ y=0.072m,z=0.034m y=0.072m,z=0.034m ） 地下厂房洞室群的地下厂房洞室群的三维围岩稳定分析三维围岩稳定分析 FEPGFEP

50、G系统在水电工程中的应用系统在水电工程中的应用工程背景工程背景进行索风营电站的建设和进一步优化索风营电站地下工程的设计和施工。通过模拟计算揭示出索风营围岩稳定状态，进行围岩稳定评价，为设计的优化和施工程序的安排提供指导。 地下厂房结构图地下厂房结构图 模拟计算的主要困难模拟计算的主要困难1.考虑横观各向同性的弹塑性；2.主变室及尾水隧洞所处岩层大部分为薄层；3.大量支护锚杆的存在如何处理；4.岩体渗流对围岩稳定的影响；岩石的横观各向同性特征岩石的横观各向同性特征 利用FEPG生成横观各向同性弹塑性计算程序软弱结构面三维造型软弱结构面三维造型 利用组合网格技术处理软弱面加锚杆岩体单元加锚杆岩体单

51、元 利用等效方法构造加锚杆的岩体单元渗流计算网格图渗流计算网格图 利用生（死）节点法处理渗流自由面 Step15_厂房位移云图厂房位移云图 Step15_厂房变形示意图厂房变形示意图Step15_厂房第一主应力厂房第一主应力 Step15_厂房屈服区厂房屈服区 Step15_衬砌支护上的变形衬砌支护上的变形Step15_衬砌支护上的第一主应力衬砌支护上的第一主应力七、结论七、结论采用元件化程序设计方法，大大降低了有限元程序的复杂性。采用脚本语言形式的有限元语言，免去了大量繁琐的，费力耗时的有限元编程。由有限元公式自动产生有限元代码，不仅节省了编程时间，并且大大提高了程序的可读性和再用性，为有限元程序的维护和发展创造了前所未有的前景。FEPG完全基于有限元方法的基本原理（虚位移原理），不受专业领域的限制，各种有限元问题和有限元方法均可采用。FEPG是一个开放源码的，任何人都可参与开发的自由软件。联系我们联系我们http:/谢谢大家！