有限元自動生成及其應(yīng)用演示

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

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

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

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

5、一種有層次的腳本語言。有限元語言應(yīng)包括微分方程表達式，形函數(shù)表達式，算法表達式等有限元公式的書寫形式。有限元語言應(yīng)與有限元的專業(yè)語言盡可能一致，以便廣大科技工作者易于理解與接受。有限元語言要適用于各種領(lǐng)域和各種有限元方法與有限元問題。關(guān)鍵問題是如何描述有限元問題，采用什么樣的方式準確表達有限元問題及其求解算法。2、人工智能技術(shù)13用FEPG生成求解二維線彈性平面應(yīng)變問題的有限元程序六、程序生成實例14問題的方程描述平衡方程幾何方程本構(gòu)方程15虛功方程虛功方程弱形式將本構(gòu)關(guān)系代入上式即可得以位移為基本未知量的虛功方程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系統(tǒng)在復(fù)合材料中的應(yīng)用25三維彈性本構(gòu)關(guān)系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系統(tǒng)在水力學(xué)中的應(yīng)用35問題描述 如右圖所示，水在1/4環(huán)形水渠內(nèi)流動，進口處x方向單寬流量為4, y方向單寬流量為0,水深為2，出口處水面高度為4。求水渠內(nèi)的單寬流量和水邊高度。 36該問題的偏微分方程水流連續(xù)性方程：水流動量方程：37該問題的偏微分方程懸移質(zhì)不平衡輸移方程 懸移質(zhì)引起的河床變形方程 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結(jié)果圖58線性正弦穩(wěn)態(tài)渦流問題（TEAM Workshop 問題7） FEPG系統(tǒng)在電磁學(xué)中的應(yīng)用59問題描述 TEAM Workshop 問題7是用來檢驗三維線性正弦渦流場分析方法計算精度的實驗?zāi)Ｐ?/p>

21、，其結(jié)構(gòu)如右圖所示。該模型由帶有方孔的鋁板和載流線圈組成，線圈中通有隨時間按正弦規(guī)律變化的電流，其幅值為2742安匝。鋁板的電導(dǎo)率為 求渦流和磁通密度分布。 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