1、1一、通用有限元軟件的局限性 通用軟件并不通用，只能求解很少的一部分有限元問(wèn)題。需要有一支訓(xùn)練有素、組織嚴(yán)密、長(zhǎng)期穩(wěn)定的有限元軟件隊(duì)伍，不斷維護(hù)和改進(jìn)其通用有限元程序。只能求解今天已有的微分方程，不能求解明天的有限元問(wèn)題,只能采用今天已出現(xiàn)的算法，不可能采用明天才有的算法。不公開(kāi)源代碼，對(duì)用戶(hù)來(lái)說(shuō)程序是不可改動(dòng)的，極大地限制了用戶(hù)的創(chuàng)造性。2二、自己編程的種種困難有限元程序復(fù)雜，代碼龐大。編程者缺乏軟件專(zhuān)業(yè)訓(xùn)練，代碼可讀性差，再用性差，程序難以維護(hù)，程序缺乏嚴(yán)密組織與發(fā)展。大多數(shù)工程師、教授與研究人員，他們主要是研究有限元方法或把有限元方法應(yīng)用于工程設(shè)計(jì)與科學(xué)研究，而不是致力于有限元軟件的開(kāi)發(fā)

2、，因此他們不可能花很多時(shí)間編寫(xiě)程序。3三、開(kāi)發(fā)FEPG的動(dòng)機(jī)與目的許多專(zhuān)家希望有自己的有限元程序，但不愿花許多時(shí)間編程。FEPG應(yīng)該不受領(lǐng)域的限制，各種專(zhuān)業(yè)的有限元問(wèn)題都可采用本系統(tǒng)求解。本系統(tǒng)不可能象通用有限元軟件那樣，通過(guò)不斷增加單元模塊和算法模塊來(lái)逐步擴(kuò)大其求解范圍。不管什么樣的有限元問(wèn)題，不管采用什么樣的算法，都能采用本系統(tǒng)求解。本系統(tǒng)是一種公開(kāi)源代碼的自由軟件，任何非本系統(tǒng)開(kāi)發(fā)者都可維護(hù)和修改。4四、由有限元公式自動(dòng)產(chǎn)生有限元程序有限元是求解偏微分方程的數(shù)值方法。有限元公式包括： （1）弱形式的微分方程表達(dá)式。 （2）形函數(shù)表達(dá)式。 （3）算法表達(dá)式。 a. 形成代數(shù)方程組的矩陣表達(dá)

3、式。 b. 計(jì)算步驟與迭代控制。5五、軟件實(shí)現(xiàn)思想1、元件化程序設(shè)計(jì)方法 把有限元程序分解為五個(gè)元件程序。6數(shù)據(jù)轉(zhuǎn)換程序可視化程序后處理程序前處理程序后處理計(jì)算程序單元計(jì)算程序單元子程序代數(shù)方程組求解程序有限元計(jì)算部分調(diào)用有限元程序結(jié)構(gòu)7有限元程序流程框圖線性穩(wěn)態(tài)問(wèn)題流程前處理程序START元件程序E元件程序SOLV求解器U元件程序DEFER元件程序計(jì)算結(jié)果圖形顯示8有限元程序流程框圖 線性動(dòng)態(tài)問(wèn)題流程前處理程序START元件程序BFT元件程序E元件程序SOLV求解器U元件程序DEFER元件程序時(shí)間計(jì)算是否結(jié)束計(jì)算結(jié)果圖形顯示否 是9有限元程序流程框圖非線性穩(wěn)態(tài)問(wèn)題流程前處理程序START元件

4、程序E元件程序SOLV求解器U元件程序迭代是否結(jié)束計(jì)算結(jié)果圖形顯示否 是DEFER元件程序10有限元程序流程框圖非線性動(dòng)態(tài)問(wèn)題流程前處理程序START元件程序BFT元件程序E元件程序SOLV求解器U元件程序DEFER元件程序迭代是否結(jié)束計(jì)算結(jié)果圖形顯示否 是時(shí)間計(jì)算是否結(jié)束 是否11自動(dòng)生成系統(tǒng)由近百個(gè)元件程序組成。元件程序是以文件名為輸入輸出參數(shù)的可獨(dú)立運(yùn)行的程序。由腳本（批命令）文件把元件程序組成一個(gè)“程序系統(tǒng)”。大大降低程序的復(fù)雜性，大大提高程序的可讀性和再用性。實(shí)現(xiàn)軟件自由化（即自由軟件）的一種軟件設(shè)計(jì)方法。12提供一種用以描述有限元問(wèn)題的語(yǔ)言，以下稱(chēng)為有限元（描述）語(yǔ)言。有限元語(yǔ)言是

5、一種有層次的腳本語(yǔ)言。有限元語(yǔ)言應(yīng)包括微分方程表達(dá)式，形函數(shù)表達(dá)式，算法表達(dá)式等有限元公式的書(shū)寫(xiě)形式。有限元語(yǔ)言應(yīng)與有限元的專(zhuān)業(yè)語(yǔ)言盡可能一致，以便廣大科技工作者易于理解與接受。有限元語(yǔ)言要適用于各種領(lǐng)域和各種有限元方法與有限元問(wèn)題。關(guān)鍵問(wèn)題是如何描述有限元問(wèn)題，采用什么樣的方式準(zhǔn)確表達(dá)有限元問(wèn)題及其求解算法。2、人工智能技術(shù)13用FEPG生成求解二維線彈性平面應(yīng)變問(wèn)題的有限元程序六、程序生成實(shí)例14問(wèn)題的方程描述平衡方程幾何方程本構(gòu)方程15虛功方程虛功方程弱形式將本構(gòu)關(guān)系代入上式即可得以位移為基本未知量的虛功方程16需填寫(xiě)的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三維正交各向異性材料計(jì)算

9、FEPG系統(tǒng)在復(fù)合材料中的應(yīng)用25三維彈性本構(gòu)關(guān)系26其中27其中28需填寫(xiě)的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平面二維水沙運(yùn)動(dòng)問(wèn)題FEPG系統(tǒng)在水力學(xué)中的應(yīng)用35問(wèn)題描述 如右圖所示，水在1/4環(huán)形水渠內(nèi)流動(dòng)，進(jìn)口處x方向單寬流量為4, y方向單寬流量為0,水深為2，出口處水面高度為4。求水渠內(nèi)的單寬流量和水邊高度。 36該問(wèn)題的偏微分方程水流連續(xù)性方程：水流動(dòng)量方程：37該問(wèn)題的偏微分方程懸移質(zhì)不平衡輸移方程 懸移質(zhì)引起的河床變形方程 38算子分裂法求解該問(wèn)題采用算子分裂法分五步驟求解上述平面二維水沙運(yùn)動(dòng)方程 1解水流擴(kuò)散方程和水位 39需填寫(xiě)的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算子分裂法求解該問(wèn)題2用基于牛頓法的最小二乘法求解水流對(duì)流方程 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算子分裂法求解該問(wèn)題3解挾沙方程擴(kuò)散項(xiàng) 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算子分裂法求解該問(wèn)題4解挾沙方程對(duì)流

19、項(xiàng) 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算子分裂法求解該問(wèn)題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結(jié)果圖58線性正弦穩(wěn)態(tài)渦流問(wèn)題（TEAM Workshop 問(wèn)題7） FEPG系統(tǒng)在電磁學(xué)中的應(yīng)用59問(wèn)題描述 TEAM Workshop 問(wèn)題7是用來(lái)檢驗(yàn)三維線性正弦渦流場(chǎng)分析方法計(jì)算精度的實(shí)驗(yàn)?zāi)Ｐ?/p>

21、，其結(jié)構(gòu)如右圖所示。該模型由帶有方孔的鋁板和載流線圈組成，線圈中通有隨時(shí)間按正弦規(guī)律變化的電流，其幅值為2742安匝。鋁板的電導(dǎo)率為 求渦流和磁通密度分布。 60求解方程及邊界條件61需填寫(xiě)的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計(jì)算渦流的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計(jì)算磁通密度的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