第三章 线弹性静力学问题

1、第三章 线弹性静力学问题第一节 用向导AppWizard解线弹性静力学问题本例将求解一个简支的深梁上表面受均布载荷作用下的平面应变问题，求出梁的变形和应力分布。图3.1.1 受均布压力的简支深梁该问题的偏微分方程表达式及定解条件如下：1）平衡方程 2）几何方程 3）本构方程 4）边界条件第二类边界条件具体操作步骤如下：设定工作路径。点击File菜单下的WorkDir选项，在弹出的对话框中填入工作路径。如图3.1.2所示图3.1.2 设定路径对话框点击AppWizard菜单，弹出公式库列表，系统提供了数学、固体力学、岩土、热传导、渗流、流体力学等领域的各种问题的公式。如图3.1.3所示。点击so

2、lid，在Project Name中填入“S”；点击“Next”进入下一步。图3.1.3 选择学科领域选择坐标系。点击左边的“2dxy”，选中二维直角坐标系，如图3.1.4所示。再点击“Next”进入下一步。图3.1.4 选择坐标系选择问题类型，即该问题是线弹性还是弹塑性。点击“elastic”，选择线弹性，如图3.1.5所示。点击“Next”进入下一步。图3.1.5 选择问题类型选择方程类型，即方程是平面应力问题还是平面应变问题，是静态问题还是动态问题；对于动态问题还可选择Newmark格式求解。点击“Static_Plane_Strain_Problem”，选择选择静态平面应力问题的方程，

3、如图3.1.6所示。点击“Next”进入下一步。图3.1.6 选择方程类型选择单元类型。如图3.1.7所示。点中“q4”，再点击中间的“Add”按钮，加入到右边方框表示采用四节点矩形单元。点击“Next”进入下一步。图3.1.7 选择单元类型选择边界单元。对于本算例，因为是应力边界条件，所以需要边界单元。点中“ll2”，再点击中间的“Add”按钮，加入到右边方框表示采用两节点线单元。如图3.1.8所示。点击“Next”进入下一步。图3.1.8 选择边界单元选择求解器。系统提供有：对称求解器（sin）、非对称求解器（nin）、Gauss-Seidel迭代求解器（gs）、逐次超松弛迭代求解器（so

4、r）、共轭梯度法求解器（cgm）、不完全LU分解求解器（ilu）等。点击“sin”，选择对称求解器。如图3.1.9所示。点击“Next”进入下一步。图3.1.9 选择求解器选择存储方式。对于有限元程序自动生成系统网络版（IFEPG）仅提供用外存的存储方式。对FEPG3.0提供外存和内存两种存储方式。选择计算数据。该数据给出求解区域的有限元数据，如图3.1.10所示。用户应注意：点击左边文件列表选择数据时，看一下右边方框中说明的方程类型和单元类型与你前面选择的是否一致，不一致系统将给出提示，一致时才可以正确运行。点击“sdispe_q4”，点击“Next”进入下一步。图3.1.10 选择有限元计

5、算数据自动生成全部有限元程序。点击“Run”即可生成全部有限元程序。如图3.1.11和图3.1.12所示。图3.1.11 点击“Run”自动生成图3.1.12 全部有限元程序计算结果插入res01.bmp图3.1.13 算例的有限元计算结果第二节 改变施加的载荷 生成全部有限元计算程序之后，在FEPG界面的左边是WorkSpace列表，它将所有生成的文件以及生成这些文件所需的原始文件都列在其中。FEPG施加荷载的方式有三种，一是分布荷载，作为材料属性给出，二是面力，通过边界单元材料属性给出，三是节点力，通过节点给出。改变施加的载荷需修改前处理中的DAT文件，以下分别按三种加载分别叙述具体如何修

6、改。a. 改变分布力的大小在WorkSpace列表中找到“preprocessor”文件夹，在修改文件DAT文件之前先看一看MTI文件，该文件定义了有限元计算中需要的所有数据表格。单击“preprocessor”文件夹中的s.mti文件，在界面的右边则显示出该文件。找到mate段如下maten pe pv fx fyi4 E10.4 E10.4 E10.4 E10.4mate 定义为材料参数表格名；n表示材料类型号；pe 表示弹性模量；pv 表示poisson比；fx 表示x方向的分布力；fy 表示y方向的分布力；在了解了mate定义的材料参数意义后，再看DAT文件。单击“preprocess

7、or”文件夹中的s.dat文件编辑该文件，找到mate信息段内容如下mate1;1.0e10;0.25;0.;0.;对应与s.mti文件，我们可以知道1 表示第一种材料，如果有多种材料可在下一行写“2”及相应的材料参数；1. 0e10 表示材料的弹性模量的值；0.25 表示材料的泊松比的值；第一个0. 表示x方向分布力的大小；第二个0. 表示y方向分布力的大小；弄清各个数值的意义我们就可以对它作相应的修改了，比如假设在x方向有1.0e3大小的力则mate信息段改成如下形式mate1;1.0e10;0.25;1.0e3;0.;下面看一下改变了载后的计算结果。首先在s.dat上单机鼠标右键弹出如下

8、命令mti s dat该mti的两个输入参数分别是s.mti 和 s.dat两个文件，输出是前处理程序prg.for和前处理显示程序prt.for。运行上述命令重新生成prg.for和prt.for。然后运行run_batch_file文件夹中的s.bat重新的计算结果如下插入res02.bmpb. 改变面力同样，我们先看一下s.mti文件。在s.mti文件中找到bf段如下bfn fx fyi4 E10.4 E10.4bf 定义为面力的参数表格名n 表示面力类型fx 面的切向方向力的分量的值fy 面的法向方向力的分量的值（参考局部标架说明）再看s.dat中对应的行，找到bf信息段如下bf1;0

9、.0;-1.0e6;1 表示第一种类型的面力0.0 表示切向方向上的力为零-1.0e6 表示法向方向上的力为1.0e6，方向和法向相反。下面给出两种不同的面力，一是将-1.0e6改为1.0e6，这表示算例中的梁是受向上的拉力，同上面所述保存s.dat文件并用mti命令重新生成prg.for和prt.for，再运行s.bat重新进行计算，得计算结果如下插入res03.bmp二是将bf信息段改为如下形式bf1;1.0e6;0.0;即在切向方向加上剪力，方向方向不加力。同样重新生成prg.for并进行计算。计算结果如下插入res04.bmpc. 施加节点力节点力的修改同样需在前处理中给出，让我们看一

10、看“preprocessor”中的s.mti文件和s.dat文件。先看s.mti中id段和disp段id fn idu idvi6 i6 i6id 表示该段定义约束信息n 表示节点编号idu x方向的约束信息存入idu数组idv y方向的约束信息存入idv数组disp fn u vi6 e10.4 e10.4disp 表示该段定义给定位移或荷载信息n 表示节点编号u x方向的给定位移或荷载信息存入u数组v y方向的给定位移或荷载信息存入v数组再看s.dat文件中id信息段和disp信息段id&nn=(nx+1)*(ny+1); i;1;1;i=1,nn1; -1;-1;nx+1; 1;

11、-1;&n=nn;&nn=(nx+1)*(ny+1) 定义总节点数i;1;1;I=1,nn 以配对的方括号表示循环，即表示先把所有节点的约束赋值为1。FEPG系统规定1表示节点自由，-1表示节点被约束。1; -1;-1; 表示1节点在x和y两个方向都被约束住。nx+1;1;-1; 表示nx+1节点仅在y方向被约束住，x方向是自由的。&n=nn 表示把最大节点后赋值给n，用户应注意加上此行，因为这将使生成的prg.for中把约束信息写入文件的循环的上限为最大节点号，也就是保证将所有的节点约束信息写入文件。再看disp信息段dispi;0.0;0.0;i=1,(nx+1)*

12、(ny+1)disp 表示该段定义给定位移和荷载信息i;0.0;0.0;i=1,(nx+1)*(ny+1) 以配对的方括号表示循环，将所有节点在x和y两个方向的位移或荷载信息都赋零。FEPG系统作如下判断：如果某节点的某个方向的约束信息是-1，即该节点在该方向被约束住，则在disp中给出的相应的值表示给定位移；如果某节点的某个方向的约束信息是1，即该节点在该方向是自由的，则在disp中给出的相应的值表示给定荷载。为了更清楚地看到施加节点荷载的结果，先将分布力和面力各方向上的荷载值置零，修改后的mate信息段和bf信息段如下mate1;1.0e10;0.25;0.0;0.;bf1;0.0;0.0

13、;在右上角节点沿水平方向向右施加大小为1.0e4的集中力，修改后的disp信息段如下dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*(ny+1);1.0e4;0.0;用mti命令重新生成前处理程序再计算后的如下结果插入res05.bmp在右上角节点沿垂直方向向上施加大小为1.0e4的集中力，修改后的disp信息段如下dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*(ny+1); 0.0;1.0e4;用mti命令重新生成前处理程序再计算后的如下结果插入res06.bmp在左上角节点沿水平方向向左施加大小为1.0e4的集中力，修改后的disp

14、信息段如下dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*ny+1; -1.0e4; 0.0;&n=(nx+1)*(ny+1);大家应注意此修改增加了一行赋值，即把总节点数赋给n，理由同在id中阐述的理由一样。为何上两例没有给这样的赋值呢？因为上两例都是给的最大节点号，所以n自然取得总节点数，否则也应有同样的赋值行。荷载正负号的规定是和总体坐标一致为正，反之为负。用mti命令重新生成前处理程序再计算后的如下结果插入res07.bmp, u1放大50000被的变形图第三节 改变约束方式上面我们已经提到id信息段给出约束信息，下面我们就举几个例子说明如何修改约束

15、信息假设物体左端固定，右端悬空，且在右上角受一大小为1.0e4向下的集中力，如下图所示。对s.dat文件中id信息段做如下修改id&nn=(nx+1)*(ny+1);i;1;1;i=1,nn(i-1)*(nx+1)+1;-1;-1;i=1,ny+1&n=nn;相应的节点荷载也作响应修改，修改后内容如下dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*(ny+1);0.0;-1.0e4;用mti命令重新生成前处理程序再计算后的如下结果插入res08.bmp，变形放大5000倍的结果假设物体左端固定，右端有一竖向支撑，且在右上角受一大小为1.0e4向下的集

16、中力，如下图所示。在上例的基础上只要增加右下角节点在y方向上的约束即可，修改后的id信息段如下id&nn=(nx+1)*(ny+1);i;1;1;i=1,nn(i-1)*(nx+1)+1;-1;-1;i=1,ny+1nx+1;1;-1;&n=nn;用mti命令重新生成前处理程序再计算后的如下结果插入res09.bmp，变形放大50000倍的结果最后一例假设物体两端固定，且在中间位置受一大小为1.0e4向下的集中力，如下图所示。修改id信息段如下id&nn=(nx+1)*(ny+1);i;1;1;i=1,nn(i-1)*(nx+1)+1;-1;-1;i=1,ny+1i*(

17、nx+1);-1;-1;i=1,ny+1&n=nn;disp信息段也作如下修改dispi;0.0;0.0;i=1,(nx+1)*(ny+1)(nx+1)*ny+(nx+2)/2;0.0;-1.0e4;&n=(nx+1)*(ny+1);用mti命令重新生成前处理程序再计算后的如下结果插入res10.bmp，变形放大50000倍的结果第四节 改变材料参数前面在修改分布力是已经提到mate信息段，并解释了各参数的意义，其中前两个参数是线弹性材料的弹性模量个泊松比，下面就修改材料参数，看计算结果如何变化。在两段固定，中间受以集中力的情况下，看不同的材料的变形情况。下面分别给出在弹性模量

18、E分别为1.0e9，5.0e9和1.0e10时mate信息段修改后的内容mate1;1.0e9;0.25;0.0;0.;mate1;5.0e9;0.25;0.0;0.;mate1;1.0e10;0.25;0.0;0.;每次修改后保存s.dat文件，用mti命令重新生成前处理程序再计算后的得如下计算结果的比较插入res11.bmp，变形放大5000倍的结果第五节 两种不同的材料III下面考虑两种不同材料时，如何用FEPG求解，假设求解下列两层材料深梁在集中力作用下的变形，如下图所示。如图中的两层材料弹性模量，泊松比，厚度相同，右上角受垂直向下大小为1.0e4大小的力，分析在此情景下深梁的变形。为

19、用FEPG解算该问题，需在前处理中作相应的修改。看一看“preprocessor”文件夹中s.mti文件的q4段，内容如下q4n nod1 nod2 nod3 nod4 matei6 i6 i6 i6 i6 i6q4 表示四节点矩形单元n 表示单元编号nod1 表示单元第一节点编号存入的数组名nod2 表示单元第二节点编号存入的数组名nod3 表示单元第三节点编号存入的数组名nod4 表示单元第四节点编号存入的数组名mate 表示该单元对应的材料号对于本算例需修改s.dat文件实现分层划分单元，上层mate的值给1，下层mate的值给2，mate对应的信息段分别给出两中材料的材料参数，修改后的

20、具体内容如下。对q4信息段作如下修改q4nx*(j-1)+i;(j-1)*(nx+1)+i;(j- 1)*(nx+1)+i+1;j*(nx+1)+i+1;j*(nx+1)+i;1;i=1,nx;j=1,ny/2nx*(j-1)+i;(j-1)*(nx+1)+i;(j- 1)*(nx+1)+i+1;j*(nx+1)+i+1;j*(nx+1)+i;2;i=1,nx;j=ny/2+1,nymate信息段也需相应修改mate1;1.0e10;0.25;0.0;0.;2;2.0e10;0.30;0.0;0.;用mti命令重新生成前处理程序再计算后得如下结果insert res12.bmp and res

21、13.bmp 变形放大5000倍的结果第六节 线弹性力学方程的PDE文件FEPG通过用户填写的PDE文件生成单元子程序。在了解如何填写PDE文件之前先回顾一下线弹性力学所对应的偏微分方程。线弹性力学偏微分方程包括平衡方程、几何方程和本构方程。在二维直接坐标系下，线弹性平面应变问题的方程具体可写成如下形式平衡方程： （1）几何方程 （2）本构方程 （3）对(1)第一个方程乘，第二个方程乘，将两式对体积积分并相加得 （4）对(4)式分布积分可得 （5）在考虑下列Cauchy公式 （6）其中和分别表示应力边界上力矢量在x和y方向的分量。将Cauchy公式和几何方程（2）式代入（5）式得 （7）将本构

22、方程代入（7）式得最终的虚功方程的弱形式 （8）其中 为一常系数有了(8)式我们就可以开始填写PDE文件，先考虑（8）式中的体积积分项，由此构造体单元。至于(8)中的边界积分项我们在下一节给出详细说明。按照（8）式给出如下的PDE文件，参阅“sdispe”文件夹中的sdispe.pde文件右边给出相应的注解disp u,v, 给出未知函数名U，Vcoor x,y, 给出总体坐标系下的坐标变量名X，Yfunc = funa,funb,func, 给出需要用到的函数名Shap %1 %2 给出单元形状类型符和节点个数gaus %3 给出每个方向上积分点的个数或单元形状类型符mass %1 给出单元

23、形状类型符load = fu fv 给出载荷向量，对应于$c6 pe = prmt(1) 从PRMT数组中取出弹性模量，$c6表示插入FORTRAN源程序$c6 pv = prmt(2) 从PRMT数组中取出泊松比$c6 fu = prmt(3) 从PRMT数组中取出X方向体力$c6 fv = prmt(4) 从PRMT数组中取出Y方向体力$c6 fact = pe/(1.+pv)/(1.-2.*pv) 由弹模和泊松比组合成系数FACT Func 给出一些函数的定义funa=+u/x 定义FUNA为，u/x表示 funb=+v/y 定义FUNB为，v/y表示 func=+u/y+v/x 定义F

24、UNC为stif 给出刚度矩阵dist =+funa;funa*fact*(1.-pv) funa;funa的第一个FUNA+funa;funb*fact*(pv) 表示，第二个FUNA表示+funb;funa*fact*(pv) ；表示积分。+funb;funb*fact*(1.-pv) +func;func*fact*(0.5-pv)load =+u*fu+v*fv 关于PDE文件还需作如下说明a. 单元形状类型符，q表示矩形单元；b. guas信息段后给的参数是数字表示每个方向上积分点的个数，如后跟单元形状类型符则表示采用节点积分；c. prmt( )数组保存了我们前面所说的前处理数据中

25、的材料参数，在PDE文件里通过赋值取出来。在了解了PDE文件如何得来以及PDE文件各内容的意义后，让我们看一看如何由PDE文件生成单元子程序，又如何连接到单元计算的E程序中。将sdispe.pde文件中的shap、gaus、mass信息段进行以下修改shap q 4gaus qmass q其它内容不变。上述修改表示将生成四节点矩形单元子程序，体积分采用节点积分。保存sdispe.pde文件，单击FEPG界面主菜单中的“Run”，弹出FEPG命令对话框，在输入栏输入pdeges sdispe seuq4其中pdeges为系统命令，其功能是由PDE文件生成GES文件和单元子程序，它的第一个参数为P

26、DE文件主名，第二个参数为将要生成的GES文件也是单元子程序文件的主名。第二个参数也可不写，此时表示生成的GES文件或单元子程序文件和PDE文件有相同的主名。另PDE文件和GES文件的关系是：将PDE文件中的shap段和gaus段用形函数库和高斯积分库作替换即得GES文件。上述命令表示由sdispe.pde文件生成seuq4.ges和seuq4.for。输入上述命令后回车，FEPG即可生成用户所需的单元子程序。把新生成的单元子程序连接到单元计算的E程序中，使用FEPG的let命令，具体做法是在“sdispe”文件夹上单击鼠标右键弹出let命令，选择并单击该命令即可实现把单元子程序连接到单元计算

27、的E程序中。第七节 边界单元在上面讲解线弹性力学方程PDE文件时，我们得到了线弹性力学虚功方程的弱形式。上面所说的PDE文件的给出了虚功方程弱形式中体积积分的描述，这里我们将考虑如何描述边界积分项。FEPG是通过添加边界单元来描述边界积分的，边界单元的作用相当于加力的边界条件或弹簧边界，这就等效于偏微分方程的第二类或第三类边界。下面我们看边界单元的PDE文件。disp u,v, 给出未知函数名u，vcoor x, 给出单元局部坐标系下的坐标变量名Shap l 2 给出单元形状类型符和节点个数，l 表示线单元， 2 表示两节点单元gaus 2 给出每个方向上积分点的个数$c6 fx=prmt(1

28、) 从前处理中取出切向表面力， $c6表示插入FORTRAN源程序$c6 fy=prmt(2) 从前处理中取出法向表面力Stif 给出刚度矩阵dist=u;u*0.0 此算例为力的边界，所以这里刚度矩阵等于零load=+u*fx+v*fy 给出载荷向量，对应于end 结束符同样使用pdeges命令由PDE文件生成相应的单元子程序。假设上面的文件名为loadbnd.pde，单击FEPG界面主菜单中的“Run”，弹出FEPG命令对话框，在输入栏输入pdeges loadbnd seull2按回车，系统即可由loadbnd.pde文件生成seull2.ges文件和边界单元的单元计算的FORTRAN源

29、程序文件seull2.for。把新生成的单元子程序连接到单元计算的E程序中，同样使用FEPG的let命令，具体做法是在“sdispe”文件夹上单击鼠标右键弹出let命令，选择并单击该命令即可实现把单元子程序连接到单元计算的E程序中。下面我们再更进一步考虑有弹簧边界的情景，即有等效于偏微分方程的第三类边界条件的边界。此时边界力矢量和和边界上的位移有关，假设 （9）其中，分别为切向和方向的弹簧刚度。将和的表达式代入（8）式得有第三类边界的弱虚功方程（10）不难看出（10）式比（8）多一项，此即由第三类边界引起的刚度矩阵的修正。下面我们再给出有第三类边界的边界单元的PDE文件的填写disp u,v,

30、 给出未知函数名u，vcoor x, 给出单元局部坐标系下的坐标变量名Shap l 2 给出单元形状类型符和节点个数，l 表示线单元， 2 表示两节点单元gaus 2 给出每个方向上积分点的个数$c6 fx=prmt(1) 从前处理中取出切向表面力， $c6表示插入FORTRAN源程序$c6 fy=prmt(2) 从前处理中取出法向表面力$c6 ax=prmt(3) 从前处理中取出弹簧切向刚度$c6 ay=prmt(4) 从前处理中取出弹簧法向刚度Stif 给出刚度矩阵dist=u;u*ax+v;v*ay 给出刚度矩阵的修正，对应于load=+u*fx+v*fy 给出载荷向量，对应于end 结

31、束符用FEPG的pdeges命令就可由此文件生成有弹簧边界的边界单元子程序，再用let命令将单元子程序连接到E程序。这里还有一点需要说明的是需在前处理中多定义两个参数，具体做法是：a. 修改前处理中s.mti文件bf段，修改后的内容如下bfn fx fy ax ayi4 E10.4 E10.4 E10.4 E10.4b. 修改前处理中s.dat文件bf段，修改后的内容如下bf1;0.0;0.0;1.0e8;2.0e8; 后两个给出弹簧切向刚度和法向刚度下图显示第5节中算例在加入弹簧边界的结果插入res14.bmp 放大5000后的变形图第八节 坐标变换边界单元定义的未知函数是局部坐标下的未知函

32、数，而最后求解的整体坐标下的未知函数，所以需把局部坐标下的未知函数转换到整体坐标下再求解。假设在边界单元中定义的局部坐标下的位移为，整体坐标下的位移是，则用下式将局部坐标下的位移转变到整体坐标下的位移 （11）其中为坐标转换矩阵。在程序中是如何实现这种坐标变换的呢？FEPG在用户添好边界单元的PDE文件并生成相应的单元子程序之后，通过getglt命令即可自动生成坐标变换的glt文件。下面对命令的使用加以说明getglt命令的格式如下getglt ls gs dim T第一个参数ls表示单元子程序文件名(ges文件主名)第二个参数gs表示坐标变换文件名（glt文件主名）第三个参数dim表示转换到

33、的整体坐标系的维数，添入数字第四个参数T表示作坐标变换，不添该参数表示不作坐标变换，此情况出现在未知函数为标量的情况。对于第7节中的边界单元可以用下列命令生成坐标变换子程序getglt seull2 seugl2 2 T此命令生成seugl2.glt和seugl2.for文件，其中seugl2.glt文件内容如下defigsub = seugl2lsub = seull2gdim = 2ldim = 1gvar gu1,gv1,gu2,gv2,lvar lu1,lv1,lu2,lv2,node = 2vartlu1 = T1 gu1 gv1 lu2 = T1 gu2 gv2 lv1 = T2

34、gu1 gv1 lv2 = T2 gu2 gv2然后使用FEPG的glt命令由glt文件生成相应的坐标变换子程序。具体做法是在glt文件上单击鼠标右键，弹出并使用glt命令，或者单击主菜单“Run”弹出FEPG命令列表对话框，在输入栏输入glt seugl2按回车后即可生成相应的坐标变换子程序。同样需要使用let命令将它连接到单元计算的E程序中。第九节 连接的接口IO文件在前面我们已多次提到用let命令把单元子程序连接到单元计算的E程序中，这里我们介绍在实现连接中一个很重要的文件IO文件。我们结合在第1节给出的算例中计算位移文件夹sdispe中的sdispe.io文件来说明该文件的意义。sdi

35、spe.io文件的内容如下，右边文字给出相应的注解。y表示有材料性质表格2 表示有限元程序采用两种单元seuq4 给出第一种单元的子程序名，对应于subroutine seuq4()seuq4 给出第一种单元的模块名，对应于seuq4.objseugl2 给出第二种单元的子程序名，对应于subroutine seugl2()seugl2,seull2 给出第二种单元的模块名，对应于seugl2.obj和seull2.obj上面已经说过，第二种单元之所以连接两个模块是因为边界单元是在局部坐标下定义的，需用一个模块实现向整体坐标的转换。此时我们看一看E程序esdispe.for的调用情况，esdi

36、spe.for含有如下行 goto (1,2), ityp1 call seuq4(r,coef,prmt,Es,Em,Ec,Ef,ne) goto 32 call seugl2(r,coef,prmt,Es,Em,Ec,Ef,ne) goto 33 continue可以看出该E程序中调用了seuq4()和seugl2()子程序，而在含有seugl2()子程序中又有 call seull2(y,coefr,prmt,els,elm,eld,ell,num)即对seull2()的调用，为实现这样的调用，在生成可执行文件时需连接seuq4.obj、seugl2.obj和seull2.obj三个目标

37、文件。下面由第六节中的PDE文件sdispe.pde，用下列命令生成另一个单元子程序。单击Run主菜单，在命令输入栏输入pdeges sdispe seuq4my这样由sdispe.pde生成了seuq4my.ges和seuq4my.for，并将seuq4my.for由FORTRAN编译器编译成sueq4my.obj。在sdispe文件夹上单击鼠标右键并选择Add files to folder命令，加入seuq4my.ges和seuq4my.for可以看到其内容和sueq4一样，只是子程序名和for文件主名不一样。下面修改sdispe.io，连接新生成的单元子程序。修改后的sdispe.io

38、文件如下y2seuq4myseuq4myseugl2seugl2,seull2在sdispe.nfe文件上单击鼠标右键，弹出并使用gnfepg命令重新生成E程序和U程序。生成后再次查看E程序esdispe.for中单元子程序的调用发生变化。 goto (1,2), ityp1 call seuq4my(r,coef,prmt,Es,Em,Ec,Ef,ne) goto 32 call seugl2(r,coef,prmt,Es,Em,Ec,Ef,ne) goto 33 continue所以用户可以通过修改io文件实现对子程序的灵活调用。第十节 计算应力的PDE文件应力计算的单元子程序同样也是由相

39、应的PDE文件生成的。在讲应力计算单元的PDE文件如何填写之前先介绍一下应力的计算方法。大家知道，由位移为基本未知量得到的位移解在全域是连续的，应力和应变在单元内部是连续的，而单元间是不连续的，即在单元边界上发生突跳。因此，对同一个节点，围绕它的不同的单元计算得到的应力的值是不同的。FEPG采用最小二乘算法对由位移算得的应力进行局部磨平。构造下列泛函 （12）其中为磨平后的应力，表示由已知位移计算得到的单元应力。 （13）其中表示弹性矩阵，表示已知位移由几何方程求得的应变。于是，用最小二乘算法对应力进行局部磨平等价于求（12）定义的泛函的极值。即求使 （14）而(14)式又等价于，即有 （15

40、）将（15）式改写为 (16)至此我们可以填写PDE文件对(16)进行求解，计算应力的PDE文件如下（右边是注解）coef u v 位移作为已知函数传入应力计算单元coor x y 定义坐标变量名disp sa,sb,sc, 定义待求的应力的三个分量名shap %1 %2 给出单元类型符和节点数gaus %3 给出每个方向上积分点的个数或单元形状类型符mass %1 给出单元形状类型符，若系数不为1，需在第二个参数给出load = fsa fsb fsc 给出应力函数名$c6 pe = prmt(1) 给出弹性模量，以下参数意义同位移计算的PDE文件$c6 pv = prmt(2)$c6 fx

41、 = prmt(3)$c6 fy = prmt(4)$c6 fact = pe/(1.+pv)/(1.-2.*pv)stif$cv funa=+u/x 计算x方向正应变，*/*表示对已知函数的导数$cv funb=+v/y 计算y方向正应变$c6 fsa=+funa*(1.-pv)+funb*(pv) $c6 fsa=fsa*fact 计算x方向的正应力$c6 fsb=+funa*(pv)+funb*(1.-pv)$c6 fsb=fsb*fact 计算y方向的正应力$cv func=+u/y+v/x 计算剪应变$c6 fsc=func*(0.5-pv)$c6 fsc=fsc*fact 计算剪应

42、力dist = sa;sa*0.0 系数矩阵为对角阵，在mass信息段给出，所以此处dist为零11计算主应力如何使用FEPG输出用户定义的物理量，例如主应力、按材料强度理论的折算应力等，下面我们通过一个平面应力问题说明如何输出并显示用户定义的物理量。先用AppWizard生成一个平面应力的例子，然后再做些修改，加上主应力的计算，输出并显示计算结果。具体的操作步骤是：a. 设置新路径；单击File主菜单，再单击Work Dir子菜单，弹出设置路径对话框输入新路径或用浏览选择一个已有路径；b. 进入FEPG向导；单击AppWizard主菜单，进入向导对话框；c. 选择固体力学；在对话框左边选择S

43、olid，在Project Name栏输入工程名比如“S”，单击Next进入下一步；d. 单击2dxy，选择二维直角坐标系，单击Next进入下一步；e. 单击elastic，选择线弹性，单击Next进入下一步；f. 单击Static_Plane_Stress_Problem，选择静态平面应力问题，单击Next进入下一步；g. 单击q4，再单击Add按钮，选择四节点矩形单元，单击Next进入下一步；h. 单击l2，再单击Add按钮，选择两节点线单元，单击Next进入下一步；i. 单击sin，选择对称用内存的求解器，单击Next进入下一步；j. 单击outcore，选择存储单元刚度，单击Next进

44、入下一步；k. 单击sdispd.q4，选择好算例数据，单击Next进入下一步；l. 单击Run，生成全部有限元程序；m. 单击File主菜单，单击Open Workspace，单击文件夹中的wsp文件，显示生成的程序。至此，我们点鼠标的形式生成了，计算静态线弹性平面应力问题的全部有限元程序。下面介绍如何修改以加入主应力。先看一下平面应力问题计算应力的PDE文件，右边文字给出说明。最好与平面应变问题计算应力的PDE对照阅读以加深理解。coef u v 位移作为已知函数传入应力计算单元coor x y 定义坐标变量名disp sa,sb,sc, 定义待求的应力的三个分量名shap %1 %2 给

45、出单元类型符和节点数gaus %3 给出每个方向上积分点的个数或单元形状类型符mass %1 给出单元形状类型符，若系数不为1，需在第二个参数给出load = fsa fsb fsc 给出应力函数名$c6 pe = prmt(1) 给出弹性模量$c6 pv = prmt(2) 给出泊松比$c6 fu = prmt(3) 给出x方向上体积力$c6 fv = prmt(4) 给出y方向上体积力$c6 fact = pe/(1.+pv)/(1.-pv) 给出弹性矩阵中的一个系数，只是为了写刚度信息 时方便，注意与平面应变的差别stif$cv funa=+u/x 计算x方向正应变，*/*表示对已知函数的导数$cv funb=+v/y 计算y方向正应变$c6 fsa=+fu