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

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

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

3、1）弱形式的微分方程表達式。 （2）形函數(shù)表達式。 （3）算法表達式。 a. 形成代數(shù)方程組的矩陣表達式。 b. 計算步驟與迭代控制。五、軟件實現(xiàn)思想五、軟件實現(xiàn)思想1、元件化程序設(shè)計方法、元件化程序設(shè)計方法 把有限元程序分解為五個元件程序。數(shù)據(jù)轉(zhuǎn)換程序可視化程序后 處 理程序前處理程序后處理計算程序單元計算程序單元子程序代數(shù)方程組求解程序有限元計算部分調(diào)用有限元程序結(jié)構(gòu)有限元程序結(jié)構(gòu)有限元程序流程框圖有限元程序流程框圖線性穩(wěn)態(tài)問題流程前處理程序START元件程序E元件程序SOLV求解器U元件程序DEFER元件程序計算結(jié)果圖形顯示有限元程序流程框圖有限元程序流程框圖 線性動態(tài)問題流程前處理程序

4、START元件程序BFT元件程序E元件程序SOLV求解器U元件程序DEFER元件程序時間計算是否結(jié)束計算結(jié)果圖形顯示否 是有限元程序流程框圖有限元程序流程框圖非線性穩(wěn)態(tài)問題流程前處理程序START元件程序E元件程序SOLV求解器U元件程序迭代是否結(jié)束計算結(jié)果圖形顯示否 是DEFER元件程序有限元程序流程框圖有限元程序流程框圖非線性動態(tài)問題流程前處理程序START元件程序BFT元件程序E元件程序SOLV求解器U元件程序DEFER元件程序迭代是否結(jié)束計算結(jié)果圖形顯示否 是時間計算是否結(jié)束 是否自動生成系統(tǒng)由近百個元件程序組成。元件程序是以文件名為輸入輸出參數(shù)的可獨立運行的程序。由腳本（批命令）文件

5、把元件程序組成一個“程序系統(tǒng)”。大大降低程序的復雜性，大大提高程序的可讀性和再用性。實現(xiàn)軟件自由化（即自由軟件）的一種軟件設(shè)計方法。提供一種用以描述有限元問題的語言，以下稱為有限元（描述）語言。有限元語言是一種有層次的腳本語言。有限元語言應(yīng)包括微分方程表達式，形函數(shù)表達式，算法表達式等有限元公式的書寫形式。有限元語言應(yīng)與有限元的專業(yè)語言盡可能一致，以便廣大科技工作者易于理解與接受。有限元語言要適用于各種領(lǐng)域和各種有限元方法與有限元問題。關(guān)鍵問題是如何描述有限元問題，采用什么樣的方式準確表達有限元問題及其求解算法。2、人工智能技術(shù)、人工智能技術(shù)用用FEPG生成求解二維線彈性生成求解二維線彈性平面

6、應(yīng)變問題的有限元程序平面應(yīng)變問題的有限元程序六、程序生成實例六、程序生成實例問題的方程描述平衡方程00yyyxyxxyxxfyxfyxxvyuyvxuxyyyxx,xyyyxxxyyyxxE5 . 0000101)21)(1 (幾何方程本構(gòu)方程虛功方程虛功方程虛功方程弱形式虛功方程弱形式 )()()(dvTuTdVvfufdVyxVyxVxyxyyyyyxxxx0)()( VyyyxyxxyxxdVvfyxufyx將本構(gòu)關(guān)系代入上式即可得以位移為基本未知量的虛功方程需填寫的需填寫的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 p

8、v fu fv 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/

9、xstifstifdist=+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

10、* *(0.5-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 fac

11、t = pe/(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

12、/y$cv fs_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系統(tǒng)在復合材料中的應(yīng)用三維彈性本構(gòu)關(guān)系三維彈性本構(gòu)關(guān)系xyxzyzzyxxyxyxzyzzyxDTxyTDTD其中其中2

13、222220000000000000010020002000nmmnmnmnnmmnmnmnnmTsin,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

14、文件文件dispdispstrsstrs#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 v1

15、3 v23 g23 g13 g12 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)$

16、c6 c2=cos(2.0$c6 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文件文件mat

17、rix te 6 6matrix 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

18、0(-v12/e1) 1/e2 (-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_

19、ja ge_i=+le_j* *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系統(tǒng)在水力學中的應(yīng)用問題描述問題描述 如右圖所示，水在1/4環(huán)形水渠內(nèi)流動，進口處x方向單寬流量為4, y方向單寬流量為0,水深為2，出口處水面高度為4。求水渠內(nèi)的單寬流量和水邊高度。 該問題的偏微分方程該問題

20、的偏微分方程0)()(yHVxHUtZxfHVUUgnxZgHyHUyxHUxyHUVxHUUtHU3/ 1222)()()()()(yfHVUVgnyZgHyHVyxHVxyHVVxHVUtHV3/1222)()()()()(水流連續(xù)性方程：水流連續(xù)性方程：水流動量方程：水流動量方程：該問題的偏微分方程該問題的偏微分方程懸移質(zhì)不平衡輸移方程懸移質(zhì)不平衡輸移方程 懸移質(zhì)引起的河床變形方程懸移質(zhì)引起的河床變形方程 yHSDyxHSDxtZyHSVxHSUtHStt)()()()()(*SStZbsmgHVUkS322*算子分裂法求解該問題算子分裂法求解該問題采用算子分裂法分五步驟求解上述平面二維

21、水沙運動方程 1解水流擴散方程和水位解水流擴散方程和水位 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 tran

22、scall transif exist stop 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

23、exist stop goto 1if not 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 %3mas

24、s %1 volmass %1 volmate 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=

25、dsqrt(hun* *hun+hvnhun+hvn* *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* *

26、volvol* *constconst load=+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 hvn

27、1 hun hvn un1 vn1 un vn 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 fh

28、u fhu fhvvect fhu fhu fhvfuncfunc$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+h

29、u_j* *un/x_jun/x_j* *dt+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=+

30、fhu_i;fhu_i load=+fhu_iload=+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 %

31、2shap %1 %2gaus %3gaus %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$

32、c6 u2v2=dsqrt(u$c6 u2v2=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

33、 load=+hsload=+hs* *sstarsstar* *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 vec

34、t x x y vect x x y vect 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文件

35、文件disp zbdisp zbcoor x,ycoor 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

36、* *v)v)$c6 ao=rfai$c6 ao=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 結(jié)果圖結(jié)果圖線性正弦穩(wěn)態(tài)渦流問題線性正弦穩(wěn)態(tài)渦流問題（TEAM Workshop TEAM Workshop 問題問題7 7） FEPG系統(tǒng)在電磁學中的應(yīng)

37、用問題描述問題描述 TEAM WorkshopTEAM Workshop 問題7是用來檢驗三維線性正弦渦流場分析方法計算精度的實驗模型，其結(jié)構(gòu)如右圖所示。該模型由帶有方孔的鋁板和載流線圈組成，線圈中通有隨時間按正弦規(guī)律變化的電流，其幅值為2742安匝。鋁板的電導率為 求渦流和磁通密度分布。 713.526 10 ()m求解方程及邊界條件求解方程及邊界條件110011sjjJ 方程：在渦流區(qū)在非渦流區(qū)000nj 邊界條件：，在無限遠處在導體和非導體的交界面處需填寫的需填寫的FEPG公式文件公式文件GIO文件文件VDE文件文件GCN文件文件GCN文件文件defidefia ell &a e

38、ll &b str a &b str a &c str 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 %3gau

39、s %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;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 c

40、url curlx curly curlzvect curl curlx curly curlz方程的方程的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

41、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 curlzafi.vde文件文件funcfunccurlx=+Az/y-Ay/zcurlx=+Az/y-A

42、y/zcurly=+Ax/z-Az/xcurly=+Ax/z-Az/xcurlz=+Ay/x-Ax/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*

43、 *|0.0;eu1|0.0;eu1|+u/x_i;A_i+u/x_i;A_i* *|sigma;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 az

44、r 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 =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.

45、0d0;0.0d0;0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;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*

46、 *sigmasigma$cv feyi = -ayr$cv feyi = -ayr* *omegaomega* *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 bxr bxi byr byi bzr bzicoef axr axi ayr ayi azr azi ur uicoef axr axi ayr