单元子程序自动生成

单元子程序自动生成----PDE文件的填写首先我们还是以热传导问题为例来说明如何填写PDE文件。热传导实例一个正方形物体长为1,导热系数是0.044。上下绝热，左边温度维持在0，右边维持在100，求温度分布。绝热0100该问题的偏微分方程如下其中ek＝0.044(热传导系数);q=0(内热源密度)求解过程启动FEPG，利用AppWizard（或gcn系统）生成有限元计算程序运行GetPre，修改Pre文件运行PreProc，用Fepg.Gid建模,划分网格运行批处理文件（bat文件)，进行计算运行PostProc，图形化显示结果我们采用GCN系统生成程序PDE（VDE）、FBCNFEGCN、GIO－－描叙微分方程－－描叙算法，系统提供大量供选择的算法－－描叙耦合关系Gcn系统包括虚功形式参数ek,eq

对应于区域积分PDE文件dispu,

给出未知函数名Ucoorx,y,

给出总体坐标系下的坐标变量X，Y。shapq4

给出单元的形状和节点个数，q表示四边形，4表示四节点单元gaus2

给出每个方向上积分点的个数mateek

eq4.4e-2;0.0;ek,eq为材料参数名，后面为相应缺省值Ell.pdestif信息段stif

给出刚度矩阵dist=+[u/x;u/x]*ek+[u/y;u/y]*ek

第一个u/x表示

第二个u/x表示[；]表示积分,对应于load信息段load=+[u]*eq

给出载荷向量，对应于end结束符Le.gcndefiaell&STARTsinaSOLVsinaLe.gioell#elemtypeq42dxy如何生成有限元程序GIOGCN文件名PDE目的：利用填写的PDE类型的文件，系统自动生成计算单元刚度矩阵、单元质量矩阵、单元阻尼矩阵和单元荷载向量的子程序。特点：对广义位移个数、节点个数、空间坐标维数、方程阶数（四阶以内）不限制，形函数、坐标变换、任意阶虚功方程、对称非对称均可；全部符号由用户定义、易于修改和阅读；适应科学计算各种形式的变化。需填写的文件:PDE文件FBC文件需填写的文件的说明：利用有限元方法求解偏微分方程问题，要将这个方程变成其虚功方程的弱形式。PDE类型文件的填写就是基于这个弱形式。PDE文件是用来处理弱形式的体积分项。FBC文件是用来处理弱形式的边界积分项。FBC文件的填写方式与PDE文件的基本相同，只是坐标变量要比PDE文件少一维。PDE文件的结构DEFI段FUNC段STIF段MASS段DAMP段LOAD段可放入DEFI段中defi信息段的填写disp

未知函数名coor

坐标变量名coef

系数变量名func

自定义函数名mate材料参数名材料参数缺省值shap

单元形状类型符节点个数gaus

单元形状类型符或每个方向积分点个数mass单元形状类型符单元质量密度damp单元形状类型符单元阻尼系数load表达式关于单元类型符的说明q四边形单元；4，8，9节点t三角形单元；3，6节点c六面体单元；8，20，27节点w四面体或三棱柱单元；4，10四面体6，18三棱柱l线单元；2，3节点关于defi信息段的说明coef——对应于微分方程中的一个变系数，对非线性问题和耦合问题需填，和算法文件对应；func——用户定义的函数，作用是进行函数替换，使虚功方程表达式简化；func信息段的填写给出用户自定义函数的具体表达式（它们往往是未知函数及其导数的线性组合），其作用是使stif、mass、load段的填写变得简洁。填写格式：自定义函数名=±[未知函数或其导数]*表达式导数定义[u/x]；[u/x,x]；{un/x}；表达式：任何fortran表达式，fortran函数，也可以是已知函数的导数。关于func信息段的说明可以多行，但每行开头必须是+[]或-[]；每个自定义的函数之间必须空一行；不能使用自定义函数的导数。stif信息段的填写连同后面的mass、damp和load信息段，共同反应了微分方程表达式。Stif信息段的目的给出计算刚度矩阵的虚功方程表达式填写方式：dist=±[未知函数或其导数；未知函数或其导数]*表达式关于stif信息段的说明dist表示分布矩阵，刚度矩阵只允许采用分布矩阵；未知函数可以是disp中定义的函数，也可以是func中定义的函数，但后者不能用导数方式；导数是对原坐标系坐标变量的导数；[*；*]表示内积，用“；”分隔，“；”前表示虚功方程中的未知函数或其导数，“；”后表示虚位移或虚应变。“；”前后变换未知不相等。表达式可以是任何fortran表达式，fortran函数，也可以是已知函数的导数。load信息段的填写load信息段给出单元荷载即方程右端项的表达式。填写方式：load=±[未知函数或其导数]*表达式说明：规定与stif信息段完全相同，可以多行填写，但换行后必须以+[]或-[]开头。如何在单元子程序中

插入FORTRAN源程序FEPG允许defi，func或stif三个信息段插入fortran程序$c6…..：空六列加入fortran；$c0…..：不空列加入fortran；$cv…..：在插入的程序中包括已知函数对坐标的导数（{*/*}）。关于插入FORTRAN源程序的说明采用1000以上的标号，避免重复；插入的fortran程序不允许空行，因为FEPG以空行为信息段的结束标志；如何附加

FORTRAN源程序在end关键字之后空一行写关键字fort后加fortran；如何由PDE文件生成单元子程序pde

pde文件名单元子程序文件名几个例子方程类型椭圆方程波动方程抛物方程稳态热传导方程算例问题的方程描述边界条件方程虚功方程在利用有限元方法求解上面的方程时，需要对上面的高阶导数项（这里是2阶）进行降阶变成虚功方程弱形式后在计算。一般我们采用分布积分公式来进行降阶处理。虚功方程虚功方程弱形式分布积分公式虚功方程弱形式的推导虚功方程弱形式最终形式PDE文件dispu,

给出未知函数名Ucoorx,y,

给出总体坐标系下的坐标变量X，Y。shapq4

给出单元的形状和节点个数，q表示四边形，4表示四节点单元gaus2

给出每个方向上积分点的个数mateek

eq4.4e-2;0.0;ek,eq为材料参数名，后面为相应缺省值stif信息段stif

给出刚度矩阵dist=+[u/x;u/x]*ek+[u/y;u/y]*ek

第一个u/x表示

第二个u/x表示[；]表示积分,对应于load信息段load=+[u]*eq

给出载荷向量，对应于end结束符抛物问题算例问题的方程描述抛物方程边界条件虚功方程虚功方程弱形式PDE文件dispu,

给出未知函数名ucoorx,y,

给出总体坐标系下的坐标变量x，y。shapt3

给出单元的形状和节点个数，t表示三角形，3表示三节点单元gaust

表示数值积分采用顶点积分masst1.0load(-4)stifdist=+[u/x;u/x]+[u/y;u/y]end波动方程问题算例问题的方程描述波动方程边界条件虚功方程虚功方程弱形式PDE文件dispu,

给出未知函数名ucoorx,y,

给出总体坐标系下的坐标变量x，y。shap%1%2%1,%2,%3由GIO文件的规定给出相应替换gaus%3

mass%12.0damp%11.0load(-2)stifdist=+[u/x;u/x]+[u/y;u/y]end二维直角坐标下的平面应力问题算例问题的方程描述平衡方程几何方程本构方程边界条件第一类边界第二类边界虚功方程虚功方程弱形式将本构关系代入上式即可得以位移为基本未知量的虚功方程PDE文件dispu,v,给出未知函数名U，Vcoorx,y,给出总体坐标系下的坐标变量名X，Y。func=funa,funb,func,给出需要用到的函数名Shapq4给出单元的形状和节点个数，q表示四边形，4表示四节点单元gaus2给出每个方向上积分点的个数matepe

pvfufv2e80.30.00.0$c6fact=pe/(1.+pv)/(1.-2.*pv)由弹模和泊松比组合成系数FACTfunc段Func

给出一些函数的定义funa=+[u/x]定义FUNA为，[u/x]表示

funb=+[v/y]定义FUNB为，[v/y]表示

func=+[u/y]+[v/x]定义FUNC为stif信息段stif

给出刚度矩阵dist=+[funa;funa]*fact*(1.-pv)+[funa;funb]*fact*(pv)[funa;funa]的第一个FUNA+[funb;funa]*fact*(pv)表示，第二个FUNA表示+[funb;funb]*fact*(1.-pv)[；]表示积分。+[func;func]*fact*(0.5-pv)load信息段end结束符load=+[u]*fu+[v]*fv给出载荷向量，对应于边界单元FBC文件dispu,v,给出未知函数名U，Vcoorx,给出单元局部坐标系下的坐标变量名Shapl2给出单元的形状和节点个数，l表示线段，2表示二节点单元gaus2给出每个方向上积分点的个数matefx

fy0.01e3stif信息段Stif

给出刚度矩阵dist=[u;u]*0.0

没有第三类边界条件，所以这里刚度矩阵等于零load信息段end结束符load=+[u]*fx+[v]*fy

给出载荷向量，对应应力计算计算应力的PDE文件coefuvcoorxydisp

sa,sb,sc,shap%1%2gaus%3mass%11.0对应load=fsa

fsb

fsc

matepe

pvfufv2e80.30.00.0$c6fact=pe/(1.+pv)/(1.-2.*pv)stif$cv

funa=+{u/x}计算$cv

funb=+{v/y}计算$c6fsa=+funa*(1.-pv)+funb*(pv)$c6fsa=fsa*fact计算$c6fsb=+funa*(pv)+funb*(1.-pv)$c6fsb=fsb*fact计算$cv

func=+{u/y}+{v/x}计算$c6fsc=func*(0.5-pv)$c6fsc=fsc*fact计算dist=[sa;sa]*0.0stif段VDE文件的填写特点：VDE文件是PDE文件的功能扩展，增加了用户定义向量和矩阵的功能，可通过向量或矩阵分量直接表述偏微分方程，从而大大简化用户原先填写PDE文件的繁杂。这样写出来的表达式和在理论推导时直接以向量或张量表达完全一致。

利用系统命令,VDE文件可变成PDE文件。如何由VDE文件生成PDE文件vde

vde文件名几个例子泊松方程算例问题的方程描述泊松方程边界条件虚功方程弱形式VDE文件dispu

给出未知函数名ucoorx,y,

z

给出总体坐标系下的坐标变量x,y,z。shap%1%2gaus%3load(-6)vectxxyz用关键字vect定义向量x，其分量为x,y,z。

各分量之间用空格隔开。

stifdist=+[u/x_i;u/x_i]以分量形式表示end相应的PDE文件dispucoorx,y,zshap%1%2gaus%3load(-6)

stifdist=+[u/x;u/x]+[u/y;u/y]+[