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

有限元自动生成与其应用演示第一页，共90页。一、通用有限元软件的局限性

通用软件并不通用，只能求解很少的一部分有限元问题。

需要有一支训练有素、组织严密、长期稳定的有限元软件队伍，不断维护和改进其通用有限元程序。

只能求解今天已有的微分方程，不能求解明天的有限元问题,只能采用今天已出现的算法，不可能采用明天才有的算法。

不公开源代码，对用户来说程序是不可改动的，极大地限制了用户的创造性。第二页，共90页。二、自己编程的种种困难有限元程序复杂，代码庞大。编程者缺乏软件专业训练，代码可读性差，再用性差，程序难以维护，程序缺乏严密组织与发展。大多数工程师、教授与研究人员，他们主要是研究有限元方法或把有限元方法应用于工程设计与科学研究，而不是致力于有限元软件的开发，因此他们不可能花很多时间编写程序。第三页，共90页。三、开发FEPG的动机与目的许多专家希望有自己的有限元程序，但不愿花许多时间编程。FEPG应该不受领域的限制，各种专业的有限元问题都可采用本系统求解。本系统不可能象通用有限元软件那样，通过不断增加单元模块和算法模块来逐步扩大其求解范围。不管什么样的有限元问题，不管采用什么样的算法，都能采用本系统求解。本系统是一种公开源代码的自由软件，任何非本系统开发者都可维护和修改。第四页，共90页。四、由有限元公式自动产生有限元程序有限元是求解偏微分方程的数值方法。有限元公式包括：（1）弱形式的微分方程表达式。（2）形函数表达式。（3）算法表达式。a.形成代数方程组的矩阵表达式。b.计算步骤与迭代控制。第五页，共90页。五、软件实现思想1、元件化程序设计方法把有限元程序分解为五个元件程序。第六页，共90页。数据转换程序可视化程序后处理程序前处理程序后处理计算程序单元计算程序单元子程序代数方程组求解程序有限元计算部分调用有限元程序结构第七页，共90页。有限元程序流程框图线性稳态问题流程前处理程序START元件程序E元件程序SOLV求解器U元件程序DEFER元件程序计算结果图形显示第八页，共90页。有限元程序流程框图

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

有限元语言是一种有层次的脚本语言。

有限元语言应包括微分方程表达式，形函数表达式，算法表达式等有限元公式的书写形式。

有限元语言应与有限元的专业语言尽可能一致，以便广大科技工作者易于理解与接受。

有限元语言要适用于各种领域和各种有限元方法与有限元问题。

关键问题是如何描述有限元问题，采用什么样的方式准确表达有限元问题及其求解算法。2、人工智能技术第十三页，共90页。用FEPG生成求解二维线弹性平面应变问题的有限元程序六、程序生成实例第十四页，共90页。问题的方程描述平衡方程几何方程本构方程第十五页，共90页。虚功方程虚功方程弱形式将本构关系代入上式即可得以位移为基本未知量的虚功方程第十六页，共90页。需填写的FEPG公式文件GIO文件VDE文件GCN文件第十七页，共90页。Strn.gcn文件defiaell&bstrastartsinasolvsinastressb第十八页，共90页。Strn.gio文件dispstrs#elemtypet32dxy第十九页，共90页。Disp.vde文件dispuvcoorxyfuncexeyexyshap%1%2gaus%3matepepvfufv1.0e10;0.3;0.0;-100.0;$c6fact=pe/(1.+pv)/(1.-2.*pv)vectuuvvectedexeyvectfsfufvmatrde22(1.-pv)pvpv(1.-pv)第二十页，共90页。Disp.vde文件funcex=+[u/x]ey=+[v/y]exy=+[u/y]+[v/x]stifdist=+[ed_i;ed_j]*de_i_j*fact+[exy;exy]*(0.5-pv)*factload=+[u_i]*fs_iend第二十一页，共90页。Disp.vde文件funcex=+[u/x]ey=+[v/y]exy=+[u/y]+[v/x]stifdist=+[ed_i;ed_j]*de_i_j*fact+[exy;exy]*(0.5-pv)*factload=+[u_i]*fs_iend第二十二页，共90页。Strs.vde文件dispsa,sb,sc,coefuvcoorxyshap%1%2gaus%3mass%1matepepvfufv1.0e10;0.3;0.0;-100.0;$c6fact=pe/(1.+pv)/(1.-2.*pv)vectuuvvectedexeyvectfsfsafsbvectsssasbscvectfffsafsbfscmatrde22(1.-pv)pvpv(1.-pv)第二十三页，共90页。Strs.vde文件stif$cvex=+{u/x}$cvey=+{v/y}$cvfs_i=+de_i_j*ed_j*fact$cvfunc=+{u/y}+{v/x}$c6fsc=func*(0.5-pv)*factdist=+[sa;sa]*0.0load=+[ss_i]*ff_iend第二十四页，共90页。三维正交各向异性材料计算——FEPG系统在复合材料中的应用第二十五页，共90页。三维弹性本构关系第二十六页，共90页。其中第二十七页，共90页。其中第二十八页，共90页。需填写的FEPG公式文件GIO文件VDE文件GCN文件第二十九页，共90页。GCN文件defiaell&bstrastartsinasolvsinastressb第三十页，共90页。Gio文件dispstrs#elemtypew43dxyz第三十一页，共90页。Disp.vde文件dispuvwcoorxyzfuncgexgeygezgeyzgexzgexyshap%1%2gaus%3arrayre(6,6)matee1e2e3v12v13v23g23g13g12citafxfyfz\127.6d3;10.3d3;10.3d3;0.32;0.32;0.32;6.0d3;6.0d3;6.0d3;45;0.0;0.0;5.e3;$c6pi=3.14159$c6s=sin(cita*pi/180.0)$c6c=cos(cita*pi/180.0)$c6c2=cos(2.0*cita*pi/180.0)$c6s2=sin(2.0*cita*pi/180.0)fvectle6fvectge6vectgevgexgeygezgeyzgexzgexyvectuuvwvectfufxfyfz第三十二页，共90页。Disp.vde文件matrixte66c*cs*s000s2s*sc*c000(-s2/2)000100000c(-s)0000sc0(-s2/2)s2/2000c2matrRe661/e1(-v12/e1)(-v13/e1)000(-v12/e1)1/e2(-v23/e2)000(-v13/e1)(-v23/e2)1/e30000001/g230000001/g130000001/g12第三十三页，共90页。Disp.vde文件func@ldeform.xyzflexyzuvw@age_i=+[le_j]*te_i_j@wgevge@lrevmatrin6re

stifdist=+[gev_i;gev_j]*re_i_jload=+[u_i]*fu_iend第三十四页，共90页。平面二维水沙运动问题——FEPG系统在水力学中的应用第三十五页，共90页。问题描述

如右图所示，水在1/4环形水渠内流动，进口处x方向单宽流量为4,y方向单宽流量为0,水深为2，出口处水面高度为4。求水渠内的单宽流量和水边高度。

第三十六页，共90页。该问题的偏微分方程水流连续性方程：水流动量方程：第三十七页，共90页。该问题的偏微分方程悬移质不平衡输移方程

悬移质引起的河床变形方程

第三十八页，共90页。算子分裂法求解该问题采用算子分裂法分五步骤求解上述平面二维水沙运动方程

1．解水流扩散方程和水位

第三十九页，共90页。需填写的FEPG公式文件GIO文件VDE文件GCN文件第四十页，共90页。GCN文件defiashola&bsholbcsholcdsholdeshole

startsinastartsinbstartsincstartsindstartsinecalltrans第四十一页，共90页。ifexiststopdelstop:1bftsolvsinacopyunodunodaifexistenddelend:2solvsinbifnotexistendgoto2solvsincsolvsindsolvsinecallpostifnotexiststopgoto1第四十二页，共90页。GIO文件fldflcsddsdczbh

#elemtypeq42dxy第四十三页，共90页。FLD.VDE文件disphu,hv,zcoorx,ycoefhun,hvn,zn,hn,funcdivvecthuhuhvvectxxyshap%1%2gaus%3mass%1volmateepsilonegenfxfy0.1d-19.8d02.0d-20.0d00.0d0vectflfxfy第四十四页，共90页。FLD.VDE文件func$c6vol=1.0d0$c6const=1.0div=+[hu/x]+[hv/y]

stif$cvu2v2=dsqrt(hun*hun+hvn*hvn)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*const

load=+[hu_i]*vol*fl_i

end第四十五页，共90页。算子分裂法求解该问题2．用基于牛顿法的最小二乘法求解水流对流方程

第四十六页，共90页。FLC.VDE文件disphuhvcoorxyfuncfhufhvcoefhun1hvn1hunhvnun1vn1unvnhn1materou1.0shap%1%2gaus%3mass%11.0vecthunhunhvnvectxxyvectfhun1fhun1fhvn1vectununvnvectun1un1vn1vecthuhuhvvectfhufhufhv第四十七页，共90页。func$cvfhun1=+{un_j/x_j}*hun+{hun/x_j}*un_j$cvfhun1=+fhun1*dt+hun1$cvfhvn1=+{un_j/x_j}*hvn+{hvn/x_j}*un_j$cvfhvn1=+fhvn1*dt+hvn1fhu=+[hu]+[hu_j]*{un/x_j}*dt+[hu_j/x_j]*un*dt+[hu]*{un_j/x_j}*dt+[hu/x_j]*un_j*dt

fhv=+[hv]+[hu_j]*{vn/x_j}*dt+[hu_j/x_j]*vn*dt+[hv]*{un_j/x_j}*dt+[hv/x_j]*un_j*dt

stifdist=+[fhu_i;fhu_i]

load=+[fhu_i]*fhun1_i

end第四十八页，共90页。算子分裂法求解该问题3．解挟沙方程扩散项

第四十九页，共90页。SDD.VDE文件disphscoorx,ycoefhn1,u,vfuncgradhsxgradhsyvectgradhsgradhsxgradhsyvectxxyshap%1%2gaus%3mass%1volmateedtrhorousrfaiomegekeg\0.1d-126501.01.01.5d-42.2d-19.8第五十页，共90页。func$c6vol=1.0d0gradhsx=+[hs/x]

gradhsy=+[hs/y]

stif$c6u2v2=dsqrt(u*u+v*v)$c6ao=rfai*omega$c6sstar=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第五十一页，共90页。算子分裂法求解该问题4．解挟沙方程对流项

第五十二页，共90页。SDC.VDE文件disphscoorxyfuncfhscoefhs1un1vn1hn1hnshap%1%2gaus%3materho1.0vectxxyvectun1un1vn1第五十三页，共90页。funcfhs=+[hs]+[hs/x_j]*un1_j*dt+[hs]*{un1_j/x_j}*dt

stifdist=+[fhs;fhs]

load=+[fhs]*(hs1-rho*(hn-hn1))

end第五十四页，共90页。算子分裂法求解该问题5．解河床变形方程

第五十五页，共90页。ZBH.VDE文件dispzbcoorx,ycoefhs,u,v,hshap%1%2gaus%3mass%1rhosmaterhosrfaiomegekegem\1.4d11.0d01.5d-42.2d-19.8d07.6d-1

第五十六页，共90页。stif$c6u2v2=dsqrt(u*u+v*v)$c6ao=rfai*omeg$cvsstar=ek*(u2v2**3.0d0/(eg*h*omeg))**emdist=+[zb;zb]*0.0d0

load=+[zb]*ao*(hs/h-sstar)

end第五十七页，共90页。结果图第五十八页，共90页。线性正弦稳态涡流问题

（TEAMWorkshop问题7）

——FEPG系统在电磁学中的应用第五十九页，共90页。问题描述

TEAMWorkshop问题7是用来检验三维线性正弦涡流场分析方法计算精度的实验模型，其结构如右图所示。该模型由带有方孔的铝板和载流线圈组成，线圈中通有随时间按正弦规律变化的电流，其幅值为2742安匝。铝板的电导率为求涡流和磁通密度分布。

第六十页，共90页。求解方程及边界条件第六十一页，共90页。需填写的FEPG公式文件GIO文件VDE文件GCN文件第六十二页，共90页。GCN文件defiaell&bstra&cstra&startcasolvcastressbstressc第六十三页，共90页。GIO文件afimagemagb第六十四页，共90页。方程的afi.vde文件dispAx,Ay,Az,ucoorx,y,zfunccurlx,curly,curlz,divshap%1%2gaus%3mass%11.0d0matesigmaepsilonomegafmufjxfjyfjz\0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;vectAAxAyAzvectxxyzvectfjfjxfjyfjzvectcurlcurlxcurlycurlz第六十五页，共90页。方程的afi.vde文件dispAx,Ay,Az,ucoorx,y,zfunccurlx,curly,curlz,divshap%1%2gaus%3mass%11.0d0matesigmaepsilonomegafmufjxfjyfjz\0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;vectAAxAyAzvectxxyzvectfjfjxfjyfjzvectcurlcurlxcurlycurlz第六十六页，共90页。afi.vde文件funccurlx=+[Az/y]-[Ay/z]curly=+[Ax/z]-[Az/x]curlz=+[Ay/x]-[Ax/y]div=+[Ax/x]+[Ay/y]+[Az/z]第六十七页，共90页。afi.vde文件stif$c6eu1=omega*sigma$c6eu2=1/fmudist=+[curl_i;curl_i]*|eu2;0.0|+[div;div]*|eu2;0.0|+[A_i;A_i]*|0.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|end第六十八页，共90页。计算涡流的mage.pde文件dispexrexieyreyiezrezicoefaxraxiayrayiazraziuruicoorxyzshap%1%2gaus%3mass%1load=fexrfexifeyrfeyifezrfezimatesigmaepsilonomegafmufjxfjyfjz\0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;第六十九页，共90页。mage.pde文件stif$cvfexr=+axi*omega*sigma-{ur/x}*sigma$cvfexi=-axr*omega*sigma-{ui/x}*sigma$cvfeyr=+ayi*omega*sigma-{ur/y}*sigma$cvfeyi=-ayr*omega*sigma-{ui/y}*sigma$cvfezr=+azi*omega*sigma-{ur/z}*sigma$cvfezi=-azr*omega*sigma-{ui/z}*sigmadist=+[exr;exr]*0.0end第七十页，共90页。计算磁通密度的magb.pde文件dispbxrbxibyrbyibzrbzicoefaxraxiayrayiazraziuruicoorxyzshap%1%2gaus%3mass%1load=fbxrfbxifbyrfbyifbzrfbzimatesigmaepsilonomegafmufjxfjyfjz\0.0d0;0.0d0;50.d0;1.256d-6;0.0d0;0.0d0;0.0d0;第七十一页，共90页。magb.pde文件stif$cvfbxr=+{azr/y}-{ayr/z}$cvfbxi=+{azi/y}-{ayi/z}$cvfbyr=+{axr/z}-{azr/x}$cvfbyi=+{axi/z}-{azi/x