第三章线弹性静力学问题精.docx

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

2、3.1.3所示。点击solid，在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）、逐

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

5、.1.10 选择有限元计算数据自动生成全部有限元程序。点击“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文件。

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

8、at上单机鼠标右键弹出如下命令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中对应的行，找

9、到bf信息段如下bf1;0.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

11、;nx+1; 1;-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)*(ny

12、+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*(nx+1);-1;-1;i=1,ny+1&n=nn;d

17、isp信息段也作如下修改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信息段，并解释了各参数的意义，其中前两个参数是线弹性材料的弹性模量个泊松比，下面就修改材料参数，看计算结果如何变化。在两段固定，中间受以集中力的情况下，看不同的材料的变形情况。下面分别给出在弹性模量E分别为1.0e9，5.0e9和1.0e10时mate信息段修改后的内

18、容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倍的结果第五节 两种不同的材料 下面考虑两种不同材料时，如何用FEPG求解，假设求解下列两层材料深梁在集中力作用下的变形，如下图所示。如图中的两层材料弹性模量，泊松比，厚度相同，右上角受垂直向下大小为1.0e4大小的力，分析在此情景下深梁的变形。为用FEPG解算该问题，需在前处理中作相应的修改。看一看“preproces

19、sor”文件夹中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对应的信息段分别给出两中材料的材料参数，修改后的具体内容如下。对q4信息段作如下修改q4nx*(j-1)+i;(j-1)*

20、(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 res13.bmp 变形放大5000倍的结果第六节 线弹性力学方程的PDE文件F

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

22、式我们就可以开始填写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 给出单元形状类型符load = fu fv 给出载荷向量，对应于$c6 pe =

23、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 定义FUNC为stif 给出刚度矩阵dist =+funa;funa*fact*

24、(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( )数组保存了我们前面所说的前处理数据中的材料参数，在PDE文件里通过赋值取出来。在了解了PDE文件如何得来以及P

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

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

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

28、6 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源程序文件seull2.for。把新生成的单元子程序连接到单元计算的E程

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

30、名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 结束符用FEPG的pdeges命令就可由此文件生成有弹簧边界的边界单元子程序，再用let命令将单元子程序连接到E程序。这里还有一点需要说明的是需在前处理中多定义两个参数，具体做法是：a. 修改前处理中s.mti文件bf段，修改后的内容如下bfn fx