有限元编程大作业报告
哈工大有限元大作业
作业一一．计算程序和结果展示1.程序clearsyms a b c x E lD=E*pi*(b*x+c)^4/64;B(:,:,1)=[-6/l^2+12*x/l^3 4/l-6*x/l^2];B(:,:,2)=[6/l^2-12*x/l^3 ,2/l-6*x/l^2];n=0;for i=1:2for j=1:2n=n+1;f=B(:,:,i)*D*transpose(B(:,:,j));k(:,:,n)=int(f,x,0,l);endendk11=k(:,:,1);k12=k(:,:,2);k21=k(:,:,3);k22=k(:,:,4);K=[k11 k12;k21 k22];K=simple(K);2.结果（1）bx+cK =[ (3*pi*E*(11*b^4*l^4 + 49*b^3*c*l^3 + 84*b^2*c^2*l^2 + 70*b*c^3*l +35*c^4))/(560*l^3), -(pi*E*(19*b^4*l^4 + 84*b^3*c*l^3 + 147*b^2*c^2*l^2 + 140*b*c^3*l + 105*c^4))/(1120*l^2), -(3*pi*E*(11*b^4*l^4 + 49*b^3*c*l^3 + 84*b^2*c^2*l^2 + 70*b*c^3*l + 35*c^4))/(560*l^3), -(pi*E*(47*b^4*l^4 + 210*b^3*c*l^3 + 357*b^2*c^2*l^2 + 280*b*c^3*l + 105*c^4))/(1120*l^2)][ -(pi*E*(19*b^4*l^4 + 84*b^3*c*l^3 + 147*b^2*c^2*l^2 + 140*b*c^3*l +105*c^4))/(1120*l^2), (pi*E*(3*b^4*l^4 + 14*b^3*c*l^3 + 28*b^2*c^2*l^2 +35*b*c^3*l + 35*c^4))/(560*l), (pi*E*(19*b^4*l^4 + 84*b^3*c*l^3 + 147*b^2*c^2*l^2 +140*b*c^3*l + 105*c^4))/(1120*l^2), (pi*E*(13*b^4*l^4 + 56*b^3*c*l^3 +91*b^2*c^2*l^2 + 70*b*c^3*l + 35*c^4))/(1120*l)][ -(3*pi*E*(11*b^4*l^4 + 49*b^3*c*l^3 + 84*b^2*c^2*l^2 + 70*b*c^3*l +35*c^4))/(560*l^3), (pi*E*(19*b^4*l^4 + 84*b^3*c*l^3 + 147*b^2*c^2*l^2 + 140*b*c^3*l +105*c^4))/(1120*l^2), (3*pi*E*(11*b^4*l^4 + 49*b^3*c*l^3 + 84*b^2*c^2*l^2 +70*b*c^3*l + 35*c^4))/(560*l^3), (pi*E*(47*b^4*l^4 + 210*b^3*c*l^3 + 357*b^2*c^2*l^2 + 280*b*c^3*l + 105*c^4))/(1120*l^2)][ -(pi*E*(47*b^4*l^4 + 210*b^3*c*l^3 + 357*b^2*c^2*l^2 + 280*b*c^3*l +105*c^4))/(1120*l^2), (pi*E*(13*b^4*l^4 + 56*b^3*c*l^3 + 91*b^2*c^2*l^2 +70*b*c^3*l + 35*c^4))/(1120*l), (pi*E*(47*b^4*l^4 + 210*b^3*c*l^3 + 357*b^2*c^2*l^2 +280*b*c^3*l + 105*c^4))/(1120*l^2), (pi*E*(17*b^4*l^4 + 77*b^3*c*l^3 +133*b^2*c^2*l^2 + 105*b*c^3*l + 35*c^4))/(560*l)]（2）ax^2+bx+c（将程序中D的直径换成“ax^2+bx+c”）K =[ (pi*E*(518*a^4*l^8+2233*a^3*b*l^7+2420*a^3*c*l^6+3630*a^2*b^2*l^6+7920*a^2*b*c*l^5 + 4356*a^2*c^2*l^4 +2640*a*b^3*l^5+8712*a*b^2*c*l^4+9702*a*b*c^2*l^3+3696*a*c^3*l^2 + 726*b^4*l^4 + 3234*b^3*c*l^3 + 5544*b^2*c^2*l^2 + 4620*b*c^3*l +2310*c^4))/(12320*l^3), -(pi*E*(938*a^4*l^8+4004*a^3*b*l^7+4290*a^3*c*l^6+6435*a^2*b^2*l^6+13860*a^2*b*c*l^ 5+7524*a^2*c^2*l^4+4620*a*b^3*l^5+15048*a*b^2*c*l^4+16632*a*b*c^2*l^3+6468*a*c^3*l ^2 +1254*b^4*l^4+5544*b^3*c*l^3+9702*b^2*c^2*l^2+9240*b*c^3*l+6930*c^4))/(73920*l^2), -(pi*E*(518*a^4*l^8 + 2233*a^3*b*l^7 + 2420*a^3*c*l^6 + 3630*a^2*b^2*l^6 +7920*a^2*b*c*l^5 + 4356*a^2*c^2*l^4 + 2640*a*b^3*l^5 + 8712*a*b^2*c*l^4 +9702*a*b*c^2*l^3 + 3696*a*c^3*l^2 + 726*b^4*l^4 + 3234*b^3*c*l^3 + 5544*b^2*c^2*l^2 + 4620*b*c^3*l + 2310*c^4))/(12320*l^3), -(pi*E*(2170*a^4*l^8 + 9394*a^3*b*l^7 +10230*a^3*c*l^6 + 15345*a^2*b^2*l^6 + 33660*a^2*b*c*l^5 + 18612*a^2*c^2*l^4 +11220*a*b^3*l^5 + 37224*a*b^2*c*l^4 + 41580*a*b*c^2*l^3 + 15708*a*c^3*l^2 +3102*b^4*l^4 + 13860*b^3*c*l^3 + 23562*b^2*c^2*l^2 + 18480*b*c^3*l +6930*c^4))/(73920*l^2)][ -(pi*E*(938*a^4*l^8 + 4004*a^3*b*l^7 + 4290*a^3*c*l^6 + 6435*a^2*b^2*l^6 + 13860*a^2*b*c*l^5 + 7524*a^2*c^2*l^4 + 4620*a*b^3*l^5 + 15048*a*b^2*c*l^4 +16632*a*b*c^2*l^3 + 6468*a*c^3*l^2 + 1254*b^4*l^4 + 5544*b^3*c*l^3 + 9702*b^2*c^2*l^2 + 9240*b*c^3*l + 6930*c^4))/(73920*l^2), (pi*E*(434*a^4*l^8 + 1848*a^3*b*l^7 + 1980*a^3*c*l^6 + 2970*a^2*b^2*l^6 + 6435*a^2*b*c*l^5 + 3564*a^2*c^2*l^4 +2145*a*b^3*l^5 + 7128*a*b^2*c*l^4 + 8316*a*b*c^2*l^3 + 3696*a*c^3*l^2 + 594*b^4*l^4 + 2772*b^3*c*l^3 + 5544*b^2*c^2*l^2 + 6930*b*c^3*l + 6930*c^4))/(110880*l),(pi*E*(938*a^4*l^8 + 4004*a^3*b*l^7 + 4290*a^3*c*l^6 + 6435*a^2*b^2*l^6 +13860*a^2*b*c*l^5 + 7524*a^2*c^2*l^4 + 4620*a*b^3*l^5 + 15048*a*b^2*c*l^4 +16632*a*b*c^2*l^3 + 6468*a*c^3*l^2 + 1254*b^4*l^4 + 5544*b^3*c*l^3 + 9702*b^2*c^2*l^2 + 9240*b*c^3*l + 6930*c^4))/(73920*l^2), (pi*E*(1946*a^4*l^8 + 8316*a^3*b*l^7 +8910*a^3*c*l^6 + 13365*a^2*b^2*l^6 + 28710*a^2*b*c*l^5 + 15444*a^2*c^2*l^4 +9570*a*b^3*l^5 + 30888*a*b^2*c*l^4 + 33264*a*b*c^2*l^3 + 12012*a*c^3*l^2 +2574*b^4*l^4 + 11088*b^3*c*l^3 + 18018*b^2*c^2*l^2 + 13860*b*c^3*l +6930*c^4))/(221760*l)][ -(pi*E*(518*a^4*l^8 + 2233*a^3*b*l^7 + 2420*a^3*c*l^6 +3630*a^2*b^2*l^6 + 7920*a^2*b*c*l^5 + 4356*a^2*c^2*l^4 + 2640*a*b^3*l^5 +8712*a*b^2*c*l^4 + 9702*a*b*c^2*l^3 + 3696*a*c^3*l^2 + 726*b^4*l^4 + 3234*b^3*c*l^3 +5544*b^2*c^2*l^2 + 4620*b*c^3*l + 2310*c^4))/(12320*l^3), (pi*E*(938*a^4*l^8 + 4004*a^3*b*l^7 + 4290*a^3*c*l^6 + 6435*a^2*b^2*l^6 + 13860*a^2*b*c*l^5 +7524*a^2*c^2*l^4 + 4620*a*b^3*l^5 + 15048*a*b^2*c*l^4 + 16632*a*b*c^2*l^3 +6468*a*c^3*l^2 + 1254*b^4*l^4 + 5544*b^3*c*l^3 + 9702*b^2*c^2*l^2 + 9240*b*c^3*l + 6930*c^4))/(73920*l^2), (pi*E*(518*a^4*l^8 + 2233*a^3*b*l^7 +2420*a^3*c*l^6 + 3630*a^2*b^2*l^6 + 7920*a^2*b*c*l^5 + 4356*a^2*c^2*l^4 +2640*a*b^3*l^5 + 8712*a*b^2*c*l^4 + 9702*a*b*c^2*l^3 + 3696*a*c^3*l^2 + 726*b^4*l^4 + 3234*b^3*c*l^3 + 5544*b^2*c^2*l^2 + 4620*b*c^3*l + 2310*c^4))/(12320*l^3),(pi*E*(2170*a^4*l^8 + 9394*a^3*b*l^7 + 10230*a^3*c*l^6 + 15345*a^2*b^2*l^6 +33660*a^2*b*c*l^5 + 18612*a^2*c^2*l^4 + 11220*a*b^3*l^5 + 37224*a*b^2*c*l^4 +41580*a*b*c^2*l^3 + 15708*a*c^3*l^2 + 3102*b^4*l^4 + 13860*b^3*c*l^3 +23562*b^2*c^2*l^2 + 18480*b*c^3*l + 6930*c^4))/(73920*l^2)][ -(pi*E*(2170*a^4*l^8 + 9394*a^3*b*l^7 + 10230*a^3*c*l^6 + 15345*a^2*b^2*l^6 +33660*a^2*b*c*l^5 + 18612*a^2*c^2*l^4 + 11220*a*b^3*l^5 + 37224*a*b^2*c*l^4 +41580*a*b*c^2*l^3 + 15708*a*c^3*l^2 + 3102*b^4*l^4 + 13860*b^3*c*l^3 +23562*b^2*c^2*l^2 + 18480*b*c^3*l + 6930*c^4))/(73920*l^2), (pi*E*(1946*a^4*l^8 +8316*a^3*b*l^7 + 8910*a^3*c*l^6 + 13365*a^2*b^2*l^6 + 28710*a^2*b*c*l^5 +15444*a^2*c^2*l^4 + 9570*a*b^3*l^5 + 30888*a*b^2*c*l^4 + 33264*a*b*c^2*l^3 +12012*a*c^3*l^2 + 2574*b^4*l^4 + 11088*b^3*c*l^3 + 18018*b^2*c^2*l^2 + 13860*b*c^3*l + 6930*c^4))/(221760*l), (pi*E*(2170*a^4*l^8 + 9394*a^3*b*l^7 + 10230*a^3*c*l^6 +15345*a^2*b^2*l^6 + 33660*a^2*b*c*l^5 + 18612*a^2*c^2*l^4 + 11220*a*b^3*l^5 +37224*a*b^2*c*l^4 + 41580*a*b*c^2*l^3 + 15708*a*c^3*l^2 + 3102*b^4*l^4 +13860*b^3*c*l^3 + 23562*b^2*c^2*l^2 + 18480*b*c^3*l + 6930*c^4))/(73920*l^2),(pi*E*(2282*a^4*l^8 + 9933*a^3*b*l^7 + 10890*a^3*c*l^6 + 16335*a^2*b^2*l^6 +36135*a^2*b*c*l^5 + 20196*a^2*c^2*l^4 + 12045*a*b^3*l^5 + 40392*a*b^2*c*l^4 +45738*a*b*c^2*l^3 + 17556*a*c^3*l^2 + 3366*b^4*l^4 + 15246*b^3*c*l^3 +26334*b^2*c^2*l^2 + 20790*b*c^3*l + 6930*c^4))/(110880*l)]（3）b=0（等截面梁)K =[ (3*pi*E*c^4)/(16*l^3),-(3*pi*E*c^4)/(32*l^2),-(3*pi*E*c^4)/(16*l^3),-(3*pi*E*c^4)/(32*l^2)] [ -(3*pi*E*c^4)/(32*l^2), (pi*E*c^4)/(16*l), (3*pi*E*c^4)/(32*l^2), (pi*E*c^4)/(32*l)][ -(3*pi*E*c^4)/(16*l^3), (3*pi*E*c^4)/(32*l^2), (3*pi*E*c^4)/(16*l^3),(3*pi*E*c^4)/(32*l^2)] [ -(3*pi*E*c^4)/(32*l^2), (pi*E*c^4)/(32*l), (3*pi*E*c^4)/(32*l^2), (pi*E*c^4)/(16*l)]总结：将结果（3）与教材等截面梁刚度矩阵比较，发现表达式一样，侧面证明了程序的正确性。
(完整版)有限元大作业matlab---课程设计例子
有限元大作业程序设计学校：天津大学院系：建筑工程与力学学院专业：01级工程力学姓名：刘秀学号：\\\\\\\\\\\指导老师：连续体平面问题的有限元程序分析[题目]：如图所示的正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力，m kNp 1=，同时在沿对角线y 轴上受一对集中压力，载荷为2KN ，若取板厚1=t ，泊松比0=v 。
[分析过程]：由于连续平板的对称性，只需要取其在第一象限的四分之一部分参加分析，然后人为作出一些辅助线将平板“分割”成若干部分，再为每个部分选择分析单元。
采用将此模型化分为4个全等的直角三角型单元。
利用其对称性，四分之一部分的边界约束，载荷可等效如图所示。
[程序原理及实现]：用FORTRAN程序的实现。
由节点信息文件NODE.IN和单元信息文件ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。
模型基本信息由文件为BASIC.IN生成。
该程序的特点如下：问题类型：可用于计算弹性力学平面问题和平面应变问题单元类型：采用常应变三角形单元位移模式：用用线性位移模式载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷材料性质：弹性体由单一的均匀材料组成约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束方程求解：针对半带宽刚度方程的Gauss消元法输入文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT.IN结果文件：输出一般的结果文件DATA.OUT程序的原理如框图：（1）主要变量：ID：问题类型码，ID＝1时为平面应力问题，ID=2时为平面应变问题N_NODE：节点个数N_LOAD：节点载荷个数N_DOF：自由度，N_DOF=N_NODE*2（平面问题）N_ELE：单元个数N_BAND：矩阵半带宽N_BC：有约束的节点个数PE：弹性模量PR：泊松比PT：厚度LJK_ELE(I,3)：单元节点编号数组，LJK_ELE(I,1)，LJK_ELE(I,2)，LJK_ELE(I,3)分别放单元I的三个节点的整体编号X(N_NODE), Y(N_NODE)：节点坐标数组，X(I),Y(I)分别存放节点I的x，y 坐标值P_LJK(N_BC,3)：节点载荷数组，P_LJK(I,1)表示第I个作用有节点载荷的节点的编号，P_LJK(I,2)，P_LJK(I,3)分别为该节点沿x,y方向的节点载荷数值AK(N_DOF,N_BAND)：整体刚度矩阵AKE(6,6)：单元刚度矩阵BB(3,6)：位移……应变转换矩阵（三节点单元的几何矩阵）DD(3,3)：弹性矩阵SS(3,6)；应力矩阵RESULT_N(N_NOF)：节点载荷数组，存放节点载荷向量，解方程后该矩阵存放节点位移DISP_E(6):：单元的节点位移向量STS_ELE(N_ELE,3)：单元的应力分量STS_ND(N_NODE,3)：节点的应力分量（2）子程序说明：READ_IN：读入数据BAND_K：形成半带宽的整体刚度矩阵FORM_KE：计算单元刚度矩阵FORM_P：计算节点载荷CAL_AREA：计算单元面积DO_BC：处理边界条件CLA_DD：计算单元弹性矩阵SOLVE：计算节点位移CLA_BB：计算单元位移……应变关系矩阵CAL_STS：计算单元和节点应力（3）文件管理：源程序文件：chengxu.for程序需读入的数据文件：BASIC.IN，NODE.IN，ELEMENT.IN（需要手工生成）程序输出的数据文件：DATA.OUT（4）数据文件格式：需读入的模型基本信息文件BASIC.IN的格式如下表需读入的节点信息文件NODE.IN的格式如下表需读入的单元信息文件ELEMENT.IN的格式如下表输出结果文件DATA.OUT格式如下表[算例原始数据和程序分析]：（1）模型基本信息文件BASIC.IN的数据为1,4,6,5,31.,0.,1.1,1,0,2,1,0,4,1,1,5,0,1,6,0,11,-0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5（2）手工准备的节点信息文件NODE.IN的数据为1 0.0 2.02 0.0 1.03 1.0 1.04 0. 0.5 1.0 0.6 2.0 0.（3）手工准备的单元信息文件ELEMENT.IN的数据为1 2 3 3 0 0 0 0 1 1 1 1 0 12 4 5 5 0 0 0 0 1 1 1 1 0 25 3 2 2 0 0 0 0 1 1 1 1 0 33 5 6 6 0 0 0 0 1 1 1 1 04 （4）源程序文件chengxu.for为：PROGRAM FEM2DDIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),&RESULT_N(500),AK(500,100)D IMENSION STS_ELE(500,3),STS_ND(500,3)OPEN(4,FILE='BASIC.IN')OPEN(5,FILE='NODE.IN')OPEN(6,FILE='ELEMENT.IN')OPEN(8,FILE='DATA.OUT')OPEN(9,FILE='FOR_POST.DAT')READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOADIF(ID.EQ.1)WRITE(8,20)IF(ID.EQ.2)WRITE(8,25)20 FORMAT(/5X,'=========PLANE STRESS PROBLEM========')25 FORMAT(/5X,'=========PLANE STRAIN PROBLEM========')CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT, & IJK_ELE,X,Y,IJK_U,P_IJK)CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,& IJK_ELE,X,Y,PE,PR,PT,AK)CALL FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, & STS_ELE,STS_ND)c to putout a data fileWRITE(9,70)REAL(N_NODE),REAL(N_ELE)70 FORMAT(2f9.4)WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I),& STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE)71 FORMA T(7F9.4)WRITE(9,72)(REAL(IJK_ELE(I,1)),REAL(IJK_ELE(I,2)),&REAL(IJK_ELE(I,3)),REAL(IJK_ELE(I,3)),&STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1, N_ELE)72 FORMAT(7f9.4)cCLOSE(4)CLOSE(5)CLOSE(6)CLOSE(8)CLOSE(9)E NDcc to get the original data in order to model the problemSUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR, &PT,IJK_ELE,X,Y,IJK_U,P_IJK)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3), & P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14)REAL ND_ANSYS(N_NODE,3)READ(4,*)PE,PR,PTREAD(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC)READ(4,*)((P_IJK(I,J),J=1,3),I=1,N_LOAD)READ(5,*)((ND_ANSYS(I,J),J=1,3),I=1,N_NODE)READ(6,*)((NE_ANSYS(I,J),J=1,14),I=1,N_ELE)DO 10 I=1,N_NODEX(I)=ND_ANSYS(I,2)Y(I)=ND_ANSYS(I,3)10 CONTINUEDO 11 I=1,N_ELEDO 11 J=1,3IJK_ELE(I,J)=NE_ANSYS(I,J)11 CONTINUEN_BAND=0DO 20 IE=1,N_ELEDO 20 I=1,3DO 20 J=1,3IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J))IF(N_BAND.LT.IW)N_BAND=IW20 CONTINUEN_BAND=(N_BAND+1)*2IF(ID.EQ.1) THENELSEPE=PE/(1.0-PR*PR)PR=PR/(1.0-PR)END IFR ETURNENDcC to form the stiffness matrix of elementSUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3), & AKE(6,6), SS(6,6)CALL CAL_DD(PE,PR,DD)CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 10 I=1,3DO 10 J=1,6SS(I,J)=0.0DO 10 K=1,310 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 20 I=1,6DO 20 J=1,6AKE(I,J)=0.0DO 20 K=1,320 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PTRETURNENDcc to form banded global stiffness matrixSUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE, & PR,PT,AK)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)N_DOF=2*N_NODEDO 40 I=1,N_DOFDO 40 J=1,N_BAND40 AK(I,J)=0DO 50 IE=1,N_ELECALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DO 50 I=1,3DO 50 II=1,2IH=2*(I-1)+IIIDH=2*(IJK_ELE(IE,I)-1)+IIDO 50 J=1,3DO 50 JJ=1,2IL=2*(J-1)+JJIZL=2*(IJK_ELE(IE,J)-1)+JJIDL=IZL-IDH+1IF(IDL.LE.0) THENELSEAK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)END IF50 CONTINUERETURNENDcc to calculate the area of elementSUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=(XIJ*YIK-XIK*YIJ)/2.0RETURNENDcc to calculate the elastic matrix of elementSUBROUTINE CAL_DD(PE,PR,DD)DIMENSION DD(3,3)DO 10 I=1,3DO 10 J=1,310 DD(I,J)=0.0DD(1,1)=PE/(1.0-PR*PR)DD(1,2)=PE*PR/(1.0-PR*PR)DD(2,1)=DD(1,2)DD(2,2)=DD(1,1)DD(3,3)=PE/((1.0+PR)*2.0)RETURNENDcc to calculate the strain-displacement matrix of elementSUBROUTINE CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)DO 10 II=1,3DO 10 JJ=1,310 BB(II,JJ)=0.0BB(1,1)=Y(J)-Y(K)BB(1,3)=Y(K)-Y(I)BB(1,5)=Y(I)-Y(J)BB(2,2)=X(K)-X(J)BB(2,4)=X(I)-X(K)BB(2,6)=X(J)-X(I)BB(3,1)=BB(2,2)BB(3,2)=BB(1,1)BB(3,3)=BB(2,4)BB(3,4)=BB(1,3)BB(3,5)=BB(2,6)BB(3,6)=BB(1,5)CALL CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DO 20 I1=1,3DO 20 J1=1,620 BB(I1,J1)=BB(I1,J1)/(2.0*AE)RETURNENDcc to form the global load matrixSUBROUTINE FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3), & RESULT_N(N_DOF)DO 10 I=1,N_DOF10 RESULT_N(I)=0.0DO 20 I=1,N_LOADII=P_IJK(I,1)RESULT_N(2*II-1)=P_IJK(I,2)20 RESULT_N(2*II)=P_IJK(I,3)RETURNENDcc to deal with BC(u) (here only for fixed displacement) using "1-0" method SUBROUTINE DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100)DO 30 I=1,N_BCIR=IJK_U(I,1)DO 30 J=2,3IF(IJK_U(I,J).EQ.0)THENELSEII=2*IR+J-3AK(II,1)=1.0RESULT_N(II)=0.0DO 10 JJ=2,N_BAND10 AK(II,JJ)=0.0DO 20 JJ=2,II20 AK(II-JJ+1,JJ)=0.0END IF30 CONTINUERETURNENDcc to solve the banded FEM equation by GAUSS eliminationSUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),AK(500,100)DO 20 K=1,N_DOF-1IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1IF(N_DOF.LE.K+N_BAND-1)IM=N_DOFDO 20 I=K+1,IML=I-K+1C=AK(K,L)/AK(K,1)IW=N_BAND-L+1DO 10 J=1,IWM=J+I-K10 AK(I,J)=AK(I,J)-C*AK(K,M)20 RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K)RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)DO 40 I1=1,N_DOF-1I=N_DOF-I1IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1IF(N_BAND.LE.N_DOF-I-1)JQ=N_BANDDO 30 J=2,JQK=J+I-130 RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)40 RESULT_N(I)=RESULT_N(I)/AK(I,1)WRITE(8,50)50 FORMAT(/12X,'* * * * * RESULTS BY FEM2D * * * * *',//8X,&'--DISPLACEMENT OF NODE--'//5X,'NODE NO',8X,'X-DISP',8X,'Y-DISP') DO 60 I=1,N_NODE60 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I)70 FORMAT(8X,I5,7X,2E15.6)RETURNENDcc calculate the stress components of element and nodeSUBROUTINECAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, &STS_ELE,STS_ND)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),DD(3,3),BB(3,6), &SS(3,6),RESULT_N(N_DOF),DISP_E(6)DIMENSION STS_ELE(500,3),STS_ND(500,3)WRITE(8,10)10 FORMAT(//8X,'--STRESSES OF ELEMENT--')CALL CAL_DD(PE,PR,DD)DO 50 IE=1,N_ELECALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 20 I=1,3DO 20 J=1,6SS(I,J)=0.0DO 20 K=1,320 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 30 I=1,3DO 30 J=1,2IH=2*(I-1)+JIW=2*(IJK_ELE(IE,I)-1)+J30 DISP_E(IH)=RESULT_N(IW)STX=0STY=0TXY=0DO 40 J=1,6STX=STX+SS(1,J)*DISP_E(J)STY=STY+SS(2,J)*DISP_E(J)40 TXY=TXY+SS(3,J)*DISP_E(J)STS_ELE(IE,1)=STXSTS_ELE(IE,2)=STYSTS_ELE(IE,3)=TXY50 WRITE(8,60)IE,STX,STY,TXY60 FORMAT(1X,'ELEMENT NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)c the following part is to calculate stress components of nodeWRITE(8,55)55 FORMAT(//8X,'--STRESSES OF NODE--')DO 90 I=1,N_NODEA=0.B=0.C=0.II=0DO 70 K=1,N_ELEDO 70 J=1,3IF(IJK_ELE(K,J).EQ.I) THENII=II+1A=A+STS_ELE(K,1)B=B+STS_ELE(K,2)C=C+STS_ELE(K,3)END IF70 CONTINUESTS_ND(I,1)=A/IISTS_ND(I,2)=B/IISTS_ND(I,3)=C/IIWRITE(8,75)I,STS_ND(I,1),STS_ND(I,2),STS_ND(I,3)75 FORMAT(1X,'NODE NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)90 CONTINUERETURNENDc FEM2D programm end[算例结果]：chengxu.for所输出的数据文件DATA.OUT数据内容如下：=========PLANE STRESS PROBLEM========* * * * * RESULTS BY FEM2D * * * * *--DISPLACEMENT OF NODE--NODE NO X-DISP Y-DISP1 .000000E+00 -.525275E+012 .000000E+00 -.225275E+013 -.108791E+01 -.137363E+014 .000000E+00 .000000E+005 -.824176E+00 .000000E+006 -.182418E+01 .000000E+00--STRESSES OF ELEMENT--ELEMENT NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00ELEMENT NO.= 2STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00ELEMENT NO.= 3STX=-.108791E+01 STY=-.137363E+01 TXY= .307692E+00ELEMENT NO.= 4STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00--STRESSES OF NODE--NODE NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00NODE NO.= 2STX=-.100000E+01 STY=-.220879E+01 TXY= .249084E+00NODE NO.= 3STX=-.105861E+01 STY=-.191575E+01 TXY= .205128E+00NODE NO.= 4STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00NODE NO.= 5STX=-.970696E+00 STY=-.166667E+01 TXY= .586081E-01NODE NO.= 6STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00[结论与体会]:通过本次的课程设计，我对有限元的概念有了更加深刻的理解，同时也弥补了平时学习是疏忽的地方，充实了有限元知识。
毕业论文：有限单元法程序报告
目录程序一：平面刚架静力分析程序（PF.FOR） (17)程序二：平面三节点有限元程序 (17)程序三：四节点矩形薄板单元程序 (24)程序一：平面刚架静力分析程序（PF.FOR）已知各杆截面均为矩形，柱截面宽0.4m,高0.5m, 梁截面宽0.4m,高0.4m，各杆E=3.65×104 MPa。
图2节点、单元编号如下图3，1.2.3…..为节点号，①②③……为单元号：图3总共有13个节点，13个单元。
计算源程序如下：! PF.FOR (A program for analysis of plane frame)! Version 6.3 2004! Main program reads the control data & organizes the whole ! calculation by calling subroutines.DIMENSION W(80000)CHARACTER IDFN*20,TITLE(5)*72READ(*,'(A12)')IDFNOPEN(3,FILE=IDFN,STATUS='OLD')READ(3,'(A72)')(TITLE(M),M=1,5)READ(3,*)E,NM,NJ,NS,NLCL1=1L2=L1+NML3=L2+NML4=L3+NML11=L4+NML12=L11+NJL21=L12+NJL22=L21+NSL31=L22+NSL41=L31+6*NMCALL IOMJS(TITLE,E,NM,NJ,NS,NLC,W(L1),W(L2),W(L3),&W(L4),W(L11),W(L12),W(L21),W(L22))CALL LCVCT (NM,W(L1),W(L2),W(L31),NJ,N)CALL LCDIA (NM,N,W(L31),W(L41),W(L41),W(L41),MAXBDW,NA) L51=L41+NL52=L51+36L53=L52+NA*2L54=L53L61=L54+N*2NW=L61+6*NM-1WRITE (*,1)NA,NW1 FORMAT(/40X,'( NA=',I6,' )',/40X,'( NW=',I6,' )')CALL FORMA (E,NM,NJ,N,NA,W(L1),W(L2),W(L3),W(L4),&W(L11),W(L12),W(L31),W(L51),W(L41),W(L52)) CALL AS(NS,N,NA,W(L21),W(L41),W(L52))CALL LDLT(N,NA,W(L41),W(L52),W(L53))DO 100 LC=1,NLCREAD (3,*)NLJL62=L61+NLJL63=L62+NLJL64=L63+NLJL71=L61L81=L71+6*NMCALL B0(LC,N,NLJ,W(L54))IF(NLJ.NE.0) CALL IOLJB(N,NLJ,W(L61),W(L62),&W(L63),W(L64),W(L54))READ(3,*)NLML82=L81+NLML83=L82+NLML84=L83+NLMCALL F0(NLM,NM,W(L71))IF(NLM.NE.0) CALL IOLMFB(NM,NJ,N,NLM,W(L81),&W(L82),W(L83),W(L84),W(L1),W(L2),W(L11),&W(L12),W(L31),W(L71),W(L54))CALL BS(NS,N,W(L21),W(L22),W(L54))CALL SLVEQ(N,NA,MAXBDW,W(L41),W(L52),W(L54))CALL OJD(NJ,N,W(L54))CALL COTF(E,NM,NJ,N,W(L1),W(L2),W(L3),W(L4),&W(L11),W(L12),W(L31),W(L54),W(L71)) NW=L84+NLM-1WRITE(*,1)NA,NW100 CONTINUEWRITE(*,'(/)')ENDSUBROUTINE IOMJS(TITLE,E,NM,NJ,NS,NLC,IST,IEN,&AR,RI,X,Y,IS,VS)! Read data of members, joints, supports & print them DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),&X(NJ),Y(NJ),IS(NS),VS(NS)CHARACTER TITLE(5)*72WRITE(*,'(/)')WRITE(*,1)(TITLE(M),M=1,5)1 FORMAT(1X,A72)WRITE(*,2)E,NM,NJ,NS,NLC2 FORMAT(/13X,'The Input Data'&//10X,'The General Information'&//6X,'E',9X,'NM',5X,'NJ',5X,'NS',5X,'NLC'& /1X,1PE10.3,4I7)READ(3,*)(IST(M),IEN(M),AR(M),RI(M),M=1,NM)WRITE(*,3)3 FORMAT(/10X,'The Information of Members'&//1X,'member',2X,'start',2X,'end',9X,'A',15X,'I') WRITE(*,4)(M,IST(M),IEN(M),AR(M),RI(M),M=1,NM)4 FORMAT(1X,I4,I8,I6,1P2E16.6)READ(3,*)(X(M),Y(M),M=1,NJ)WRITE(*,5)5 FORMAT(/10X,'The Joint Coordinates'&//1X,'joint',11X,'X',17X,'Y')WRITE(*,6)(M,X(M),Y(M),M=1,NJ)6 FORMAT(1X,I4,2F18.6)READ (3,*)(IS(M),VS(M),M=1,NS)WRITE(*,7)7 FORMAT(/10X,'The Information of Supports'&//4X,'IS',9X,'VS')WRITE(*,8)(IS(M),VS(M),M=1,NS)8 FORMAT(1X,I5,F16.6)RETURNENDSUBROUTINE LCVCT(NM,IST,IEN,LV,NJ,N)! Determine location vector of elementDIMENSION IST(NM),IEN(NM),LV(6,NM)DO 100 M=1,NMI=IST(M)*3J=IEN(M)*3LV(1,M)=I-2LV(2,M)=I-1LV(3,M)=ILV(4,M)=J-2LV(5,M)=J-1LV(6,M)=J100 CONTINUEN=NJ*3RETURNENDSUBROUTINE LCDIA(NM,N,LV,MIN,IBDW,LD,MAXBDW,NA)! Determine location of diagonal elements of global stiffness ! matrix ADIMENSION LV(6,NM),MIN(N),IBDW(N),LD(N)DO 100 I=1,NMIN(I)=I100 CONTINUEDO 400 M=1,NMMINLV=LV(1,M)DO 200 I=2,6IF (LV(I,M).LT.MINLV) MINLV=LV(I,M)200 CONTINUEDO 300 I=1,6IF (MINLV.LT.MIN(LV(I,M))) MIN(LV(I,M))=MINLV300 CONTINUE400 CONTINUEMAXBDW=0DO 500 I=1,NIBDW(I)=I-MIN(I)+1IF(IBDW(I).GT.MAXBDW) MAXBDW=IBDW(I)500 CONTINUELD(1)=IBDW(1)DO 600 I=2,NLD(I)=LD(I-1)+IBDW(I)600 CONTINUENA=LD(N)RETURNENDSUBROUTINE RLCS(M,NM,NJ,IST,IEN,X,Y,RL,C,S)! Calculate length, cosine & sine of member DIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ)I=IST(M)J=IEN(M)X1=X(J)-X(I)Y1=Y(J)-Y(I)RL=SQRT(X1*X1+Y1*Y1)C=X1/RLS=Y1/RLRETURNENDSUBROUTINE KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,&X,Y,C,S,E1,E2,E3,E4)! Calculate element stiffness matrix along local axes DIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ),AR(NM),RI(NM) CALL RLCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S)E1=E*AR(M)/RLE2=12.0*E*RI(M)/(RL*RL*RL)E3=0.5*E2*RLE4=0.6666667*E3*RLRETURNENDSUBROUTINE KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)! Calculate element stiffness matrix along global axes DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),&X(NJ),Y(NJ),AE(6,6)CALL KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4) A1=E1*C*C+E2*S*SA2=(E1-E2)*C*SA3=E1*S*S+E2*C*CA4=E3*SA5=E3*CA6=E4AE(1,1)=A1AE(2,1)=A2AE(2,2)=A3AE(3,1)=-A4AE(3,2)=A5AE(3,3)=A6AE(4,1)=-A1AE(4,2)=-A2AE(4,3)=A4AE(4,4)=A1AE(5,1)=-A2AE(5,2)=-A3AE(5,3)=-A5AE(5,4)=A2AE(5,5)=A3AE(6,1)=-A4AE(6,2)=A5AE(6,3)=0.5*A6AE(6,4)=A4AE(6,5)=-A5AE(6,6)=A6RETURNENDSUBROUTINE FORMA (E,NM,NJ,N,NA,IST,IEN,AR,RI,&X,Y,LV,AE,LD,A)! Form global stiffness matrix ADIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),&LV(6,NM),AE(6,6),LD(N)DOUBLE PRECISION A(NA)DO 300 M=1,NMCALL KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)DO 200 I=1,6DO 100 J=1,IIF (LV(I,M).GE.LV(J,M)) THENA(LD(LV(I,M))-LV(I,M)+LV(J,M))&=A(LD(LV(I,M))-LV(I,M)+LV(J,M))+AE(I,J)ELSEA(LD(LV(J,M))-LV(J,M)+LV(I,M))&=A(LD(LV(J,M))-LV(J,M)+LV(I,M))+AE(I,J)END IF100 CONTINUE200 CONTINUE300 CONTINUERETURNENDSUBROUTINE AS (NS,N,NA,IS,LD,A)! Introduce support conditions into global stiffness matrix A DIMENSION IS(NS),LD(N)DOUBLE PRECISION A(NA)DO 100 M=1,NSI=3*(IS(M)/10)-3+MOD(IS(M),10)A(LD(I))=1D22100 CONTINUERETURNENDSUBROUTINE LDLT (N,NA,LD,A,T)! Solve equations (1) - decomposition of matrix A DIMENSION LD(N)DOUBLE PRECISION A(NA),T(N),SUMDO 300 I=2,NLDI=LD(I)I1=I-LDI+LD(I-1)+1DO 200 J=I1,I-1LDJ=LD(J)J1=J-LDJ+LD(J-1)+1IF(I1.GT.J1) J1=I1SUM=0.0D0DO 100 K=J1,J-1SUM=SUM+T(K)*A(LDJ-J+K)100 CONTINUET(J)=A(LDI-I+J)-SUMA(LDI-I+J)=T(J)/A(LDJ)A(LDI)=A(LDI)-T(J)*A(LDI-I+J)200 CONTINUE300 CONTINUERETURNENDSUBROUTINE SLVEQ (N,NA,MAXBDW,LD,A,B)! Solve equations (2) - forward & back substitution DIMENSION LD(N)DOUBLE PRECISION A(NA),B(N)DO 200 I=2,NLDI=LD(I)I1=I-LDI+LD(I-1)+1DO 100 J=I1,I-1B(I)=B(I)-A(LDI-I+J)*B(J)100 CONTINUE200 CONTINUEDO 300 I=1,NB(I)=B(I)/A(LD(I))300 CONTINUEDO 500 I=N-1,1,-1IMIN=I+MAXBDWIF(IMIN.GT.N) IMIN=NDO 400 J=I+1,IMINLDJ=LD(J)J1=J-LDJ+LD(J-1)+1IF(I.GE.J1) B(I)=B(I)-A(LDJ-J+I)*B(J)400 CONTINUE500 CONTINUERETURNENDSUBROUTINE B0 (LC,N,NLJ,B)! Initialize joint load vector BDOUBLE PRECISION B(N)WRITE (*,1)LC,NLJ1 FORMAT(/15X,'Loading Case',I3&//10X,'The Loadings at Joints'&//17X,'NLJ=',I4)DO 100 I=1,NB(I)=0.0D0100 CONTINUERETURNENDSUBROUTINE IOLJB (N,NLJ,ILJ,PX,PY,PM,B)! Read data of loading at joint & form joint load vector B DIMENSION ILJ(NLJ),PX(NLJ),PY(NLJ),PM(NLJ)DOUBLE PRECISION B(N)READ (3,*)(ILJ(M),PX(M),PY(M),PM(M),M=1,NLJ)WRITE (*,1)1 FORMAT(/2X,'ILJ',11X,'PX',14X,'PY',15X,'PM')WRITE (*,2)(ILJ(M),PX(M),PY(M),PM(M),M=1,NLJ)2 FORMAT(1X,I4,2F16.4,F18.5)DO 100 M=1,NLJI=ILJ(M)*3B(I-2)=PX(M)B(I-1)=PY(M)B(I)=PM(M)100 CONTINUERETURNENDSUBROUTINE F0(NLM,NM,F)! Initialize terminal forces of membersDIMENSION F(6,NM)WRITE (*,1)NLM1 FORMAT(/10X,'The Loadings at Members'&//17X,'NLM=',I4)DO 200 J=1,NMDO 100 I=1,6F(I,J)=0.0100 CONTINUE200 CONTINUERETURNENDSUBROUTINE IOLMFB(NM,NJ,N,NLM,ILM,ITL,PV,DST,&IST,IEN,X,Y,LV,F,B)! Read data of loading at member & calculate fixed-end forces, ! add equivalent joint loads to vector BDIMENSION ILM(NLM),ITL(NLM),PV(NLM),DST(NLM),IST(NM),& IEN(NM),X(NJ),Y(NJ),LV(6,NM),F(6,NM)DOUBLE PRECISION B(N)READ(3,*)(ILM(M),ITL(M),PV(M),DST(M),M=1,NLM)WRITE (*,1)1 FORMAT(/2X,'ILM',2X,'ITL',11X,'PV',12X,'DST')WRITE(*,2)(ILM(M),ITL(M),PV(M),DST(M),M=1,NLM)2 FORMAT(1X,I4,I5,F16.4,F16.6)DO 100 M=1,NLML=ILM(M)CALL RLCS (L,NM,NJ,IST,IEN,X,Y,RL,C,S) D1=DST(M)D2=RL-D1IF (ITL(M).EQ.1.OR.ITL(M).EQ.3)THENP1=PV(M)*CP2=-PV(M)*SENDIFIF(ITL(M).EQ.2.OR.ITL(M).EQ.4)THENP1=PV(M)*SP2=PV(M)*CENDIFIF(ITL(M).EQ.1.OR.ITL(M).EQ.2)THENF1=-P1*D2/RLF4=-P1-F1F2=-P2*D2*D2*(RL+2.0*D1)/(RL*RL*RL)F5=-P2-F2F3=-P2*D1*D2*D2/(RL*RL)F6=P2*D1*D1*D2/(RL*RL)ENDIFIF(ITL(M).EQ.3.OR.ITL(M).EQ.4)THENG=P2*D1*D1/(12.0*RL*RL)F3=-G*((6.0*RL-8.0*D1)*RL+3.0*D1*D1)F6=G*D1*(4.0*RL-3.0*D1)F5=-6.0*G*D1*(2.0-D1/RL)F2=-P2*D1-F5F4=-P1*D1*D1/(2.0*RL)F1=-P1*D1-F4ENDIFF(1,L)=F(1,L)+F1F(2,L)=F(2,L)+F2F(3,L)=F(3,L)+F3F(4,L)=F(4,L)+F4F(5,L)=F(5,L)+F5F(6,L)=F(6,L)+F6B(LV(1,L))=B(LV(1,L))-F1*C+F2*SB(LV(2,L))=B(LV(2,L))-F1*S-F2*CB(LV(3,L))=B(LV(3,L))-F3B(LV(4,L))=B(LV(4,L))-F4*C+F5*SB(LV(5,L))=B(LV(5,L))-F4*S-F5*CB(LV(6,L))=B(LV(6,L))-F6100 CONTINUERETURNENDSUBROUTINE BS(NS,N,IS,VS,B)! Introduce support conditions into joint load vector B DIMENSION IS(NS),VS(NS)DOUBLE PRECISION B(N)DO 100 M=1,NSI=3*(IS(M)/10)-3+MOD(IS(M),10)B(I)=VS(M)*1D22100 CONTINUERETURNENDSUBROUTINE OJD(NJ,N,B)! Print joint displacementsDOUBLE PRECISION B(N)WRITE (*,1)1 FORMAT(/13X,'The Results of Calculation'&//10X,'The Joint Displacements'&//1X,'joint',8X,'u',15X,'v',14X,'phi') WRITE (*,2)(M,B(3*M-2),B(3*M-1),B(3*M),M=1,NJ)2 FORMAT(1X,I4,1P3E16.6)RETURNENDSUBROUTINE COTF(E,NM,NJ,N,IST,IEN,AR,RI,X,Y,LV,B,F)! Calculate & print terminal forces of members DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),&LV(6,NM),F(6,NM)DOUBLE PRECISION B(N),U1,U2,U3,U4,U5,U6WRITE (*,1)1 FORMAT(/10X,'The Terminal Forces'&//1X,'member',4X,'N(st)',6X,'Q(st)',7X,'M(st)',& 6X,'N(en)',6X,'Q(en)',7X,'M(en)')DO 100 M=1,NMCALL KEBAR(M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4) U1=B(LV(1,M))*C+B(LV(2,M))*SU2=-B(LV(1,M))*S+B(LV(2,M))*CU3=B(LV(3,M))U4=B(LV(4,M))*C+B(LV(5,M))*SU5=-B(LV(4,M))*S+B(LV(5,M))*CU6=B(LV(6,M))F(1,M)=F(1,M)+E1*(U1-U4)F(2,M)=F(2,M)+E2*(U2-U5)+E3*(U3+U6)F(3,M)=F(3,M)+E3*(U2-U5)+E4*(U3+0.5*U6)F(4,M)=F(4,M)+E1*(U4-U1)F(5,M)=F(5,M)+E2*(U5-U2)-E3*(U3+U6)F(6,M)=F(6,M)+E3*(U2-U5)+E4*(0.5*U3+U6)WRITE(*,2)M,F(1,M),F(2,M),F(3,M),F(4,M),F(5,M),F(6,M) 2 FORMAT(1X,I4,2(2F11.3,F12.3))100 CONTINUERETURNEND输入数据文件************************************************** ** 114811150上机试题1 ** **************************************************3.65E7 13 13 12 11 2 0.2 4.1667E-32 3 0.2 4.1667E-33 6 0.16 2.1333E-32 5 0.16 2.1333E-34 5 0.2 4.1667E-35 6 0.2 4.1667E-36 9 0.16 2.1333E-35 8 0.16 2.1333E-38 9 0.2 4.1667E-37 8 0.2 4.1667E-38 10 0.16 2.1333E-310 12 0.16 2.1333E-311 13 0.2 4.1667E-30 00 100 2010 010 1010 2020 020 1020 2025 1030 030 1030 1011 012 013 041 042 043 071 072 073 0111 0112 0113 032 150 0 03 0 0 -50 10 0 0 -25 63 2 -250 54 4 -10 107 4 -10 108 2 -200 59 3 10 1010 3 10 10输出文件************************************************* * * * 114811150上机试题1 * * * *************************************************The Input DataThe General InformationE NM NJ NS NLC3.650E+07 13 13 12 1The Information of Membersmember start end A I1 12 2.000000E-01 4.166700E-032 23 2.000000E-01 4.166700E-033 3 6 1.600000E-01 2.133300E-034 25 1.600000E-01 2.133300E-035 4 5 2.000000E-01 4.166700E-036 5 6 2.000000E-01 4.166700E-037 6 9 1.600000E-01 2.133300E-038 5 8 1.600000E-01 2.133300E-039 8 9 2.000000E-01 4.166700E-0310 7 8 2.000000E-01 4.166700E-0311 8 10 1.600000E-01 2.133300E-0312 10 12 1.600000E-01 2.133300E-0313 11 13 2.000000E-01 4.166700E-03The Joint Coordinatesjoint X Y1 .000000 .0000002 .000000 10.0000003 .000000 20.0000004 10.000000 .0000005 10.000000 10.0000006 10.000000 20.0000007 20.000000 .0000008 20.000000 10.0000009 20.000000 20.00000010 25.000000 10.00000011 30.000000 .00000012 30.000000 10.00000013 30.000000 10.000000The Information of SupportsIS VS11 .00000012 .00000013 .00000041 .00000042 .00000043 .00000071 .00000072 .00000073 .000000111 .000000112 .000000113 .000000( NA= 258 )( NW= 927 )Loading Case 1The Loadings at JointsNLJ= 3ILJ PX PY PM2 150.0000 .0000 .000003 .0000 .0000 -50.00000 10 .0000 .0000 -25.00000The Loadings at MembersNLM= 6ILM ITL PV DST3 2 -250.0000 5.0000004 4 -10.0000 10.0000007 4 -10.0000 10.0000008 2 -200.0000 5.0000009 3 10.0000 10.00000010 3 10.0000 10.000000The Results of CalculationThe Joint Displacementsjoint u v phi1 9.727862E-21 -7.663517E-21 -5.946892E-202 8.890687E-02 -1.049797E-04 -7.120787E-033 1.470285E-01 -2.323802E-04 -7.465712E-034 9.683457E-21 -3.536778E-20 -5.930753E-205 8.886287E-02 -4.844902E-04 -7.160652E-036 1.469822E-01 -7.581843E-04 5.116445E-047 1.058868E-20 -2.196870E-20 -6.233835E-208 8.890696E-02 -3.009411E-04 -6.177452E-039 1.470137E-01 -3.792984E-04 -1.977168E-0310 8.890696E-02 -3.520154E-02 -7.782787E-0311 0.000000E+00 0.000000E+00 0.000000E+0012 8.890696E-02 -7.411547E-02 -7.782787E-0313 0.000000E+00 0.000000E+00 0.000000E+00The Terminal Forcesmember N(st) Q(st) M(st) N(en) Q(en) M(en)1 76.635 97.279 594.689 -76.635 -97.279 378.0972 93.002 -27.030 -129.904 -93.002 27.030 -140.3963 27.030 93.002 90.396 -27.030 156.998 -410.3724 25.691 -16.367 -248.192 -25.691 116.367 -415.4805 353.678 96.835 593.075 -353.678 -96.835 375.2706 199.797 45.396 110.296 -199.797 -45.396 343.6647 -18.366 42.799 66.708 18.366 57.201 -138.7178 -25.747 37.514 -70.087 25.747 162.486 -554.7759 57.201 81.634 177.624 -57.201 18.366 138.71710 219.687 155.887 706.717 -219.687 -55.887 352.15111 .000 .000 25.000 .000 .000 -25.00012 .000 .000 .000 .000 .000 .00013 .000 .000 .000 .000 .000 .000( NA= 258 )( NW= 951 )Press any key to continue程序二：平面三节点有限元程序，如图1所示的悬臂梁，受均布荷载q=1N/mm2 作用。
有限元大作业报告
有限元计算分析报告***********院&……*设计专业2008年六月目录试题一 (1)试题二 (4)试题三 (6)试题四 (8)试题一问题描述：图示为一带圆孔的单位厚度的正方形平板，在x 方向作用均布压力0.25Mpa ，试用三节点常应变单元和六节点三角形单元对平板进行有限元分析。
图1.（1）分析：（1）该平板的受力属于平面应力问题，在对称外载荷的作用下，模型的受力也是对称的。
故而对该问题也可以简化为4/1平板倒圆角的模型。
在对称均布压力作用下，平板的内部应力图也应该是对称的，平板受到沿x 方向的正应力x σ和y 方向的应力y σ，板面上没有力的作用，即：0=zy τ，0=zx τ。
垂直于板面没有正应力作用，z σ=0。
即该数学模型与z 无关，仅仅是x 、y 的函数。
（2）简化为41平板计算时，可以将倒圆左边界视作x 约束而y 方向无约束，下边界视作y 方向有约束而x 方向没有约束。
分别用三节点和六节点画出单元网格。
（3） 由于平板的中间孔存在集中应力，所以在孔的附近的有限元网格需要细化，而远离孔的网格就可以不画那么细了。
有限元网格划分结果如下图1.（2）所示；数学建模：按照1/4计算，取点（0，0，0），由矢量（0.024，0.024，0）通过surface 创建平面，再建立2DArcangle 画出90度圆弧，将平面打断删去多余部分便得到了几何模型。
六节点应力图六节点应变图三节点应力图三节点应变图三节点不同数目网格应力图三节点不同数目网格数目应变图试题二问题描述：图示 2.（1）为带方孔（边长为80mm）的悬臂梁，其上受部分均布载荷（p=10Kn/m）作用，试采用一种平面单元，对图示两种结构进行有限元分析，并就方孔的布置进行分析比较，如将方孔设计为圆孔，结果有何变化？（板厚为1mm，材料为钢）分析：（1）该题同第一个问题一样是属于平面问题，是平面的应变问题。
（2）平面及面内无z方向应力分量，且限制了z向位移。
(完整word版)有限元分析大作业报告要点
有限元分析大作业报告试题1：一、问题描述及数学建模图示无限长刚性地基上的三角形大坝，受齐顶的水压力作用，试用三节点常应变单元和六节点三角形单元对坝体进行有限元分析，并对以下几种计算方案进行比较：（1）分别采用相同单元数目的三节点常应变单元和六节点三角形单元计算；（2）分别采用不同数量的三节点常应变单元计算；（3）当选常应变三角单元时，分别采用不同划分方案计算。
该问题属于平面应变问题，大坝所受的载荷为面载荷，分布情况及方向如图所示。
二、采用相同单元数目的三节点常应变单元和六节点三角形单元计算1、有限元建模（1）设置计算类型：两者因几何条件和载荷条件均满足平面应变问题，故均取Preferences 为Structural（2）选择单元类型：三节点常应变单元选择的类型是Solid Quad 4 node182；六节点三角形单元选择的类型是Solid Quad 8 node183。
因研究的问题为平面应变问题，故对Element behavior（K3）设置为plane strain。
（3）定义材料参数：弹性模量E=2.1e11，泊松比σ=0.3（4）建几何模型：生成特征点；生成坝体截面（5）网格化分：划分网格时，拾取lineAB和lineBC，设定input NDIV 为15；拾取lineAC，设定input NDIV 为20，选择网格划分方式为Tri+Mapped，最后得到600个单元。
（6)模型施加约束：约束采用的是对底面BC 全约束。
大坝所受载荷形式为Pressure ，作用在AB 面上，分析时施加在L AB 上，方向水平向右，载荷大小沿L AB 由小到大均匀分布。
以B 为坐标原点，BA 方向为纵轴y ，则沿着y 方向的受力大小可表示为：}{*980098000)10(Y y g gh P -=-==ρρ2、 计算结果及结果分析 （1） 三节点常应变单元三节点常应变单元的位移分布图三节点常应变单元的应力分布图（2）六节点三角形单元六节点三角形单元的变形分布图六节点三角形单元的应力分布图①最大位移都发生在A点，即大坝顶端，最大应力发生在B点附近，即坝底和水的交界处，且整体应力和位移变化分布趋势相似，符合实际情况；②结果显示三节点和六节点单元分析出来的最大应力值相差较大，原因可能是B点产生了虚假应力，造成了最大应力值的不准确性。
有限元编程作业
*Elset, elset=__PickedSurf9_S4, internal, instance=Ball-1
7, 149, 162, 292, 485, 487,……,1098, 1192, 1218, 1260
*Elset, elset=__PickedSurf9_S2, internal, instance=Ball-1
*Element, type=C3D4
1, 163, 164, 165, 166
……
1357, 297, 49, 47, 51
*Nset, nset=_PickedSet2, internal, generate
1, 318, 1
*Elset, elset=_PickedSet2, internal, generate
**定义材料Mat-Ball和Mat-Plate
*Material, name=Mat-Ball
*Density
7800.,
*Elastic
2.068e+11, 0.3
*Material, name=Mat-Plate
*Density
7800.,
*Elastic
2.078e+11, 0.3
*Element, type=C3D8R
1, 243, 244, 17, 16, 1561, 1562, 1335, 1334
……
6135, 6327, 5480, 6359, 7645, 6798, 6797, 7677
**内部节点集
*Nset, nset=_PickedSet2, internal, generate
** STEP: Step-1
**定义一般静态分析步
有限元大作业试验报告
有限元大作业魏博宇力学111 3111631031一、划分单元，确定半带宽。
x13 24 6 85 7 9 11 13 15 10 12 14 1618 20 22 2417 19 21 23 25 27 29 31 33 3526 28 30 32 34 36y单元数：36 ； 节点数 28。
1.单元局部节点编号： y x ijm ijm yx12 3 4 56 7 8 9 10 11121314 1516 17 18 19 20 212223 24 25 26 27 28、单元号 1 2 3 4 5 6 7 8 9 10 11 12i 2 4 5 5 7 8 8 9 9 11 12 12j 3 5 3 6 8 5 9 6 10 12 8 13 m 1 2 2 3 4 4 5 5 6 7 7 8 单元号13 14 15 16 17 18 19 20 21 22 23 24i 13 13 14 14 16 17 17 18 18 19 19 20j 9 14 10 15 17 12 18 13 19 14 20 15 m 8 9 9 10 11 11 12 12 13 13 14 14 单元号25 26 27 28 29 30 31 32 33 34 35 36i 20 22 23 23 24 24 25 25 26 26 27 27j 21 23 17 24 18 25 19 26 20 27 21 28 m 15 16 16 17 17 18 18 19 19 20 20 212.节点坐标结点号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 x 0 0 1 0 1 2 0 1 2 3 0 1 2 3 y 0 1.5 1.5 3 3 3 4.5 4.5 4.5 4.5 6 6 6 6 结点号15 16 17 18 19 20 21 22 23 24 25 26 27 28 x 4 0 1 2 3 4 5 0 1 2 3 4 5 6 y 6 7.5 7.5 7.5 7.5 7.5 7.5 9 9 9 9 9 9 93.带状性1 5 10 15 20 25 28可求半带宽D=（7+1）× 2 = 16二、载荷向节点移置。
有限元课程大作业
金属坯料挤压过程有限元分析一、前言：金属挤压是将放在挤压模具内的金属锭坯从一端施加外力，强迫其从特定的模孔中流出，获得所需要的断面形状和尺寸的制品的技术。
冷挤压时由于材料是在冷态下成形，而且变形量一般都很大，挤压过程中作用在模具上的单位压力很大,此时模具有开裂破坏的可能，对压力机也构成威胁，金属坯料在通过模具过程中，坯料与模具之间产生相当大的应力，这就要求模具需要有相当大的强度、硬度、以及耐磨性，因此冷挤压时要进行挤压力的计算。
挤压力的计算是模具设计的重要依据，也是选择挤压设备的依据。
模具角度、接触表面的摩擦系数、坯料变形量都会影响应力变化，在保证加工要求的前提下，应当通过适当方式降低坯料及模具之间的应力。
通过有限元分析，得出应力分布图，分析变形区域、死区，对模具进行优化改进。
二、有限元介绍：ANSYS概述ANSYS软件是融结构、热、流体、电磁、声学于一体的大型通用有限元软件，可广泛地用于核工业、铁道、石油化工、航空航天、机械制造、能源、汽车交通、国防军工、电子、土木工程、生物医学、水利、日用家电等一般工业及科学研究。
该软件提供了不断改进的功能清单，具体包括：结构高度非线性分析、电磁分析、计算流体力学分析、设计优化、接触分析、自适应网格划分及利用ANSYS参数设计语言扩展宏命令功能。
ANSYS软件功能强大，主要特点有：实现多场及多场耦合分析；实现前后处理、求解及多场分析统一数据库的一体化；具有多物理场优化功能；强大的非线性分析功能;多种求解器分别适用于不同的问题及不同的硬件设备；支持异种、异构平台的网络浮动，在异种、异构平台上用户界面统一、数据文件全部兼容；强大的并行计算功能支持分布式并行及共享内存式并行；多种自动网格划分技术；良好的用户开发环境。
ANSYS不仅支持用户直接创建模型，也支持与其他CAD软件进行图形传递，其支持的图形传递有：SAT、Parasolid、STEP。
相应地，可以进行接口的常用CAD 软件有：Unigraphics、Pro/Engineer、I-Deas、Catia、CADDS、SolidEdge、SolidWorks等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
本科生实验报告书四节点等参单元有限元分析的FORTRAN程序目录1．问题概述 (1)2．四节点四边形等参单元介绍 (1)3．单元应力磨平方法介绍 (4)4．程序流程设计 (6)程序设计概述流程图5．程序结构及程序说明 (8)6．程序应用及算例分析 (9)算例概述算例ANSYS求解算例程序数值解算例分析7. 总结 (15)1. 问题概述等参单元是有限元方法中使用最广泛的单元类型。
等参单元的位移模式和坐标变换均采用相同的形函数，这种坐标变换叫做等参变换。
通过等参变换可以将自然（局部）坐标中几何形状规则的单元转换成总体（笛卡尔）坐标中形状扭曲的单元，因而使得单元有较好的适应性。
本问题首先对平面四节点四边形等参单元的形函数、应力矩阵和等效节点力矩阵、应力磨平公式等的推导和计算求解。
并通过设计FORTRAN 求解程序进行编程求解，最后给出算例（受集中荷载的悬臂梁）并进行求解，将解与ANSYS 的解进行比较。
在这个过程中，采用了高斯三点积分和高斯两点积分，这种积分方法的求解效率较高而且精度也较好。
在问题的最后，尝试去分析引起数值解误差的原因，并分析四节点等参单元的若干特性。
2. 四节点四边形等参单元介绍边长为2的正方形单元（如下图所示），在其形心处安置一个局部坐标。
单元角结点i 的坐标（,）分别为,因此单元四条边界的方程可用简单公式和逐一给出。
图2-1 母单元1 234 0图2-2 四边形单元形函数的表达式：位移函数： ∑==41i N ui i u ∑==41i N v i i v坐标变换式：∑==41i x N x i i ，∑==41i y N y i i单元应变矩阵{}[]{}[]{}ee B B B B B x v y u y v x u δδε4321==⎪⎪⎪⎭⎪⎪⎪⎬⎫⎪⎪⎪⎩⎪⎪⎪⎨⎧∂∂+∂∂∂∂∂∂=式中{}[]TTT T T 4321e δδδδδ=——单元节点的位移列阵；[]{}），，，（432100,,,,=⎭⎬⎫⎩⎨⎧=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=i v u N N N N B i i i x i y i y i x i i δ这里记号x N i x ∂∂=/N ,i ，y /N y ,i ∂∂=i N 。
根据复合函数求导规则，有[]⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧y i x i y i x i i i N N J N N y y x x N N ,,,,,,,,,,ηξηξηξ从而有[]⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧-ηξ,,1,,i i y i x i N N J N N因此，单元内的应力可以表示为{}[][][][]{}e e xy y s B D δδτσσσ==⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=x应力矩阵[]S 可以写成分块形式[][]4321S S S S S =对于平面应力情形，[][][]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---==2/)1(2/)1(1,,,,,,2x i y i yi y i x i x i i i N N N N N N E B D S μμμμμ 单元刚度矩阵是[][][][]TeAK B D B hdxdy =⎰⎰其中积分采用三点高斯积分，),(),(111111j i j nip i nipj i f d d f ηξωωζξηξ⎰⎰∑∑--===2nip n =（高斯积分点的总数），i ω和j ω是加权系数，i ξ和j η是单元内的坐标.。
对于三点高斯积分，高斯积分点的位置： 6.01-=η,110.6, 5.09.0ξω=-=，0.02=η220.0,8.0ξω==，6.03=η330.6, 5.0ξω=。
结构刚度矩阵e eK K =∑结构结点荷载列阵e eP P =∑注意，对于上两式中e∑的理解不是简单的叠加而是按照对应的自由度集成。
总刚平衡方程[]{}{}K P δ=从式上式求出：{}[]{}1K P δ-=将{}δ回代入{}[][]{}D B σδ=和{}[]{}e B δε=中，得到{}ε和{}σ。
等效节点力(1) 集中力 将集中荷载作用点取为结点，集中荷载就可以直接转化为等效结点力。
(2) 体积力 等效结点力按下式计算其中积分采用高斯2点积分，),(),(111111j i j nip i nipj i f d d f ηξωωζξηξ⎰⎰∑∑--===2nip n =（高斯积分点的总数），i ω和j ω或i W 是加权系数，i ξ和j η是单元内的坐标。
对于2点高斯积分，高斯积分点的位置：;;;(3) 表面力 单元的ij 边上两个结点的等效结点力按下式计算3. 应力磨平方法介绍为了减少改进应力结果的工作量，可以采用单元应力的局部磨平。
当单元足够小时，磨平可以在各个单元内进行。
对于等参元来说，有了单元内应力局部磨平的方法，可以方便地利用精度较高的高斯积分点(最佳应力点)的应力值来改进等参元结点应力的近似性质。
以二维单元为例，如果是二次等参元，插值函数中的完全多项式是二次，即p＝2。
对于C0型单元m＝1。
这时p - m+1=2,则在2×2个高斯积分点上有限元的应力计算值有较高的精度。
如果进行应力磨平时只要求得到4各角点的改进应力值，即使在计算位移时是8结点二次等参元，但进行应力磨平时，单元结点数可只取4个，即用二维双线性单元进行应力磨平。
磨平式各应力分量分别进行，这时应力磨平插值函数Ni'应采用双线性函数，即（）在此情况下，4个结点上的应力可由高斯点上的应力外推得到，即令在2×2高斯积分点上有。
4个高斯积分点的座标(见图3-1)如下:图3-1 二维等参单元将高斯点坐标代入式得到下面的等式：）式中等号左边的应力列阵是有限元中已求出的4个高斯点相应的应力分量；是磨平后应力的结点值；转换矩阵由式的第二式代入高斯点坐标后的插值函数值构成。
由（）式求逆可得（）其中（）各应力分量均可用（）求解。
这种改进结点应力的方法亦称之为应力插值外推。
求得改进的应力结点值后，如需要求单元内部的应力值仍可按第一式进行计算。
采用单元局部应力磨平的方法，对于同一结点，由不同相邻单元求得的应力改进值通常是不相同的。
可把相关单元求得的改进结点值再取平均作为最后的结点应力值。
4.程序流程设计4.1程序设计概述本程序的设计以P3单元的FORTRAN程序为基础，在其架构之上修改而来。
整体的架构与P3单元相似，但是在应力应变矩阵、主单元刚度矩阵、整体刚度矩阵、等效结点力以及应力磨平的算法方面有着较大的差异。
具体的算法方面，在刚度矩阵的计算用了3点高斯积分，在等效结点力的计算用了2点高斯积分。
应力磨平采用了单元内应力局部磨平的方法，对于同一结点，由不同相邻单元求得的应力改进值通常是不相同的。
把相关单元求得的改进结点值再取平均作为最后的结点应力值。
4.2流程图5.程序结构5.1程序结构5.2子程序功能说明[P4INPUT] 读入单元数据、节点数据、节点约束条件及各类单元属性[DEA] 形成主对角元素地址[P4SSM] 形成总体刚度矩阵[P4STIFB] 计算单元刚度矩阵[P4BMATB] 计算当前单元的应变矩阵[P4MODB] 计算当前单元的弹性矩阵[P4DBE] 计算应力矩阵[P4XJACM] 计算当前单元的形函数和雅可比行列式[BOUNDARY] 边界条件处理[LDLT] LDLT分解结构刚度矩阵[P4LOAK] 计算单元荷载矩阵形成节点总荷载向量用置0置1法处理边界条件时产生的荷载向量修正项，修正节点总荷载向量[SOLV] 处理荷载向量,回代求结点位移[P4STREB] 计算单元应力和主应力[P4GAUSTR] 单元内应力局部磨平[GSTREB] 计算磨平后的单元应力[P4SUMSTRS] 计算磨平后各节点的应力6.程序应用及算例分析6.1算例概述考虑一个右端固定的悬臂梁如下图5-1，其材料参数为：，，；几何参数为：L=8mm, H=2mm, b=1mm。
在左端上方有集中荷载F=-100N，不考虑悬臂梁的自重，求解悬臂梁的挠度以及各截面的应力。
图6-16.2算例ANSYS求解挠度图图6-2 X方向应力图6-3Y方向应力图6-46.3算例程序数值解将悬臂梁进行单元网格划分如下图，划分为8个单元15个结点。
图6-5将划分结果及相关参数输入程序的input文档，启动程序进行求解。
求解结果如下：表6-1 节点应力表6-2 节点位移6.4算例分析为了验证P4求解程序的准确性，我们将悬臂梁最大受拉区（顶边）的位移解和应力解在ANSYS下的和在本程序下得到的进行比较如下图：图6-6 位移解的比较由图可以看出，本程序的位移数值解与ANSYS的解十分接近，最大的节点误差在%左右。
图6-7 受拉区x向应力解比较由图可以看出，本程序的应力数值解与ANSYS的解十分接近，最大的节点误差在%左右。
应力磨平结果分析可见对于四边形等参单元应力磨平前后，局部坐标为（0，0）的点的应力值不变。
7.总结本程序采用的四结点平面等参单元具有较高的求解精度，对于平面的问题也有很好的适应性。
P4等参单元的求解还是和P3三角形常应变单元相比还是比较复杂的。
虽然求解的流程基本一致，但是等参单元的求解需要解决刚度阵和等效节点力求解时的数值积分问题，本程序中采用应用比较广泛的高斯求积法。
关于高斯积分的阶数，在李人宪所著的《有限元法基础》中指出，“用数值积分代替精确积分, 无疑计算时会产生误差. 为尽量减少这一误差,希望采用尽可能高的数值积分阶次”。
文献亦同时指出:鉴于N 很大时计算费时, 建议在二维情况下, 四节点单元的N 取为2 , 八节点单元的N 取为3 较好。
本程序则混合采用二阶和三阶高斯积分以获得较好的精确性和较高的计算效率。
至于在求解本算例和ANSYS求解存在10%左右的误差主要是因为4节点等参单元缺少边中点，不能很好地适应悬臂梁的弯曲变形，加之单元的划分比较少（总共8个单元），网格的解就存在着一定的误差。
其次，在应力磨平方面基于4个结点上的应力可由高斯点上的应力外推得到的原理，可能会存在一定的误差。