有限元编程大作业报告
哈工大有限元大作业
作业一一．计算程序和结果展示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等。
大作业报告参考2有限元学习心得
有限元学习心得吴清鸽车辆工程 50110802411短短八周的有限元课已经结束。
关于有限元，我一直停留在一个很模糊的概念。
我知道这是一个各个领域都必须涉及的点，只要有关于CAE分析的，几乎都要涉及有限元。
总体来说，这是一门非常重要又有点难度的课程。
有限元方法(finite element method) 或有限元分析(finite element analysis)，是求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要基础性原理。
将它用于在科学研究中，可成为探究物质客观规律的先进手段。
将它应用于工程技术中，可成为工程设计和分析的可靠工具。
本课程教学基本内容有固体力学和结构力学简介；有限元法基础；桁架、梁、刚架、二维固体、板和壳、三维固体的有限元法；建模技术；热传导问题的有限元分析；PATRAN软件的使用.通过有限元分析课程学习使我了解和掌握了一些有限元知识：1．简要了解二维和三维固体以及桁架、梁和板结构的三组基本力学方程，即表示位移-应变关系的几何方程，表示应力-应变关系的本构方程和表示内力-外力关系的平衡方程。
2．了解利用能量法形成有限元离散系统方程的基本原理，即哈密尔顿原理。
掌握有限元分析的基本方法及步骤，包括域的离散、位移插值、构造形函数、单元有限元方程的建立、坐标变换、整体有限元方程的组装、整体有限元方程的求解技术。
3．具体深入的了解并掌握桁架结构、梁结构、刚架结构、二维固体、板和壳结构、三维固体的有限元法分析技术，包括他们具体的形函数构造，应变矩阵，局部坐标系和整体坐标系中的单元矩阵。
各种结构的实例研究。
4．了解并掌握建立高质量建模所涉及的各种关键技术。
包括单元类型的选择，单元畸形的限制，不同阶数单元混用时网格的协调性问题，对称性的应用（平面对称、轴对称、旋转对称、重复对称），由多点约束方程形成刚域及应用（模拟偏移、不同自由度单元的连接、网格协调性的施加）等，以及多点约束方程的求解。
有限元编程大作业报告
本科生实验报告书四节点等参单元有限元分析的程序目录1．问题概述 (1)2．四节点四边形等参单元介绍 (1)3．单元应力磨平方法介绍 (4)4．程序流程设计 (6)4.1 程序设计概述4.2 流程图5．程序结构及程序说明 (8)6．程序应用及算例分析 (9)6.1 算例概述6.2 算例求解6.3 算例程序数值解6.4 算例分析7. 总结 (15)7. 附录……………………………………………………()1.问题概述等参单元是有限元方法中使用最广泛的单元类型。
等参单元的位移模式和坐标变换均采用相同的形函数，这种坐标变换叫做等参变换。
通过等参变换可以将自然（局部）坐标中几何形状规则的单元转换成总体（笛卡尔）坐标中形状扭曲的单元，因而使得单元有较好的适应性。
本问题首先对平面四节点四边形等参单元的形函数、应力矩阵和等效节点力矩阵、应力磨平公式等的推导和计算求解。
并通过设计求解程序进行编程求解，最后给出算例（受集中荷载的悬臂梁）并进行求解，将解及的解进行比较。
在这个过程中，采用了高斯三点积分和高斯两点积分，这种积分方法的求解效率较高而且精度也较好。
在问题的最后，尝试去分析引起数值解误差的原因，并分析四节点等参单元的若干特性。
2.四节点四边形等参单元介绍边长为2的正方形单元（如下图所示），在其形心处安置一个局部坐标οξη。
单元角结点i的坐标（ξi,ηi）分别为±1,因此单元四条边界的方程可用简单公式ξ±1=0和η±1=0逐一给出。
图2-1 母单元图2-2 四边形单元形函数N i的表达式：N i =(1+ξi ξ)(1+ηi η)/4位移函数：∑==41i N u i 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==⎪⎪⎪⎭⎪⎪⎪⎬⎫⎪⎪⎪⎩⎪⎪⎪⎨⎧∂∂+∂∂∂∂∂∂= 式中{}[]TT T 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 。
工作报告之ansys大作业实验报告
ansys大作业实验报告【篇一：ansys有限元分析实验报告】ansys有限元分析实验报告梅晨2013200303飞机机翼模态分析1. 问题描述对一个飞机机翼进行模态分析。
机翼沿长度方向的轮廓是一致的，横截面由直线和样条曲线定义。
机翼的一端固定在机体上，另一端悬空。
要求分析得到机翼的模态自由度。
机械的几何模型如图1所示，弹性模量取38?103pa,泊松比0.3，密度为8.3e?5kg/m3。
图 12. gui操作步骤(1)定义单元类型。
定义两种单元类型plane42和solid45，如图2所示。
图 2（2）定义材料参数。
定义ex=38000,prxy=0.3,dens=8.3e-5 如图3,图4所示。
图3图4（4）建立几何模型，首先生成关键点，然后通过关键点再生成直线。
并通过spline thru kps命令画出弧线。
最后再根据线生成机翼的截面。
各步骤如下图所示：图5图6图7图8（4）网格划分。
运用mesh tool命令对机翼进行网格划分。
各步骤如下图所示：【篇二：ansys有限元分析实验报告】ansys有限元分试验报告ansys试验报告一、 ansys简介:ansys软件是融结构、流体、电场、磁场、声场分析于一体的大型通用有限元分析软件。
由世界上最大的有限元分析软件公司之一的美国ansys开发，它能与多数cad软件接口，实现数据的共享和交换，如pro/engineer, nastran, autocad等，是现代产品设计中的高级cae工具之一。
本实验我们用的是ansys12.1软件。
二、试验题目：（6）如图所示，l/b=10，a= 0.2b ,b= (0.5-2)a，比较 b 的变化对最大应力?x 的影响；并与（5）比较。
l 三、题目分析：该问题是平板受力后的应力分析问题。
我们通过使用ansys软件求解，首先要建立上图所示的平面模型，然后在平板一段施加位移约束，另一端施加载荷，最后求解模型，用图形显示，即可得到实验结果。
有限元分析报告Ansys大作业
有限元分析作业作业名称扳手静态受力分析姓名学号宁波理工学院班级题目：扳手静态受力分析:扳手的材料参数为：弹性模量E=210GPa,泊松比u=0.3：此模型在左侧内六角施加固定位移约束，在右侧表面竖直方向上施加648 N的集中力。
10模型如下图：1-11.定义工作文件名和文件标题（1）定义工作文件名:执行File-Chang Jobname-3090601048（2）定义工作标题：执行File-Change Tile-3090601048（3）更改工作文件储存路径：执行File-Chang Directory-E:\ANSYS2.定义分析类型、单元类型及材料属性（1）定义分析类型，执行Main Menu-Preferences,如下图所示：2-1（2）定义单元类型，执行Main Menu-Preprocessor-Element Type-Add 弹出Element Type 对话框.如下图所示：2-2（3）定义材料属性执行Main menu-Preprocessor-Material Props-Material models,在Define material model behavior对话框中，双击Structual-Linear-Elastic-Isotropic.如下图所示：2-33.导入几何模型将模型导入到ANSYS，执行File-Import—PRAR…—浏览上述模型，如下图所示：3-13-24. 网格划分执行Main Menu-Preprocessor-meshing-Mesh Tool命令，考虑到零件的复杂性，采用智能网格划分，精度为1，其他选项为默认，如下图所示：4-14-25. 加载以及求解（1）添加位置约束执行Solution-apply-structural-displacement-on areas（对六角内表面进行约束），如下图所示：5-15-2（2）添加载荷，执行Solution-apply-structural-force-on keypoints，如下图所示：5-35-4（3）求解执行Main menu-Solution-Solve-Current LS,求解。
有限元程序设计大作业(DOC)
1.不同板宽的孔边应力集中问题姓名：胡宇学号：21201201282.摘要本文采用MATLAB和FOTRAN四节点平面单元，利用有限元数值解法对不同板宽的孔边应力集中问题进行了数值模拟研究。
对于不同的板宽系数，并且与解析解进行（半板宽b/孔半径r），得到了不同的应力集中系数1了比较，验证了有限元解的正确性，并且得出了解析解的适用范围。
3.引言通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、ABAQUS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。
具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。
随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中。
MATLAB是由美国MATHWORKS公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。
它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。
4.MATLAB部分1，计算模型本程序采用MATLAB编程，编制平面四边形四节点等参元程序，用以求解近似平面结构问题。
本程序的研究对象为中央开有小孔的长方形板，选取的材料参数为：板厚h=1、材料强度E=1.0e11 Pa、泊松比mu=0.3。
此外，为方便网格的划分和计算，本文所取板的长度与宽度相等。
其孔半径为r=1，板宽为2b待定。
由于本程序的目的在于验证有限元解的正确性和确定解析解的适用范围，因此要求网格足够细密，以满足程序的精度要求。
同时为了减小计算量，我采取网格径向长度递增的网格划分方法。
此种方法特点是，靠近小孔部分的网格细密，在远离小孔的过程中，网格逐渐变得稀疏。
有限元分析大作业报告要点
有限元分析大作业报告试题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）六节点三角形单元六节点三角形单元的变形分布图六节点三角形单元的应力分布图（3）计算数据表单元类型最小位移（mm）最大位移（mm）最小应力（Pa）最大应力（Pa）三节点0 0.0284 5460.7 392364六节点0 0.0292 0.001385 607043 （4）结果分析①最大位移都发生在A点，即大坝顶端，最大应力发生在B点附近，即坝底和水的交界处，且整体应力和位移变化分布趋势相似，符合实际情况；②结果显示三节点和六节点单元分析出来的最大应力值相差较大，原因可能是B点产生了虚假应力，造成了最大应力值的不准确性。
有限元大作业
有限元大作业第一篇：有限元大作业有限元应力分析报告大作业机械与运载工程学院车辆四班龙恒 20110402415 2014年8月30日一、问题描述桦木板凳材料参数如图形状参数：长40mm，宽30mm，高45mm（其他详细参数见零件图）通过施加垂直于板凳上表面的均匀载荷600N，分析板凳的应变和应力？二、使用inventor进行建模及应力分析1、通过inventor建立板凳3D模型利用草图拉伸等方法建立与零件图中尺寸一致的三维立体板凳模型2、点选环境下的应力分析开始对板凳进行应力分析3、根据所给条件设置材料等参数、将安全系数设为屈服强度，因为板凳主要受压变形点开“木材（桦木）”根据前面所给参数对其进行参数设置4、固定约束如图板凳的4个脚底面设置为固定约束，使得板凳受载后，脚底面不会沿垂直方向位移，模拟真实情况5、施加载荷在板凳上表面施加大小为600N的垂直均布载荷（这里是模拟一个成人坐上去的重力）6、划分网格通过设置网格的尺寸参数来划分出5种不同网格数量，从而得出5种不同网格数划分得出的应力应变分布图，最后分析划分不同网格数对结果的影响。
（1）网格最大（2）网格较大（3）网格一般大小（4）网格较小（5）网格最小7、求解得出结果得出5组不同网格数所得数据（应力云图，应变云图，所有结果数据）（1）网格数1437根据应力云图可知，红色地方所受的应力最大，最大应力为：15.48Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.001434μl（2）网格数8651根据应力云图可知，红色地方所受的应力最大，最大应力为：18.88Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.001755μl（3）网格数20484根据应力云图可知，红色地方所受的应力最大，最大应力为：22.62Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.002103μl（4）网格数41578根据应力云图可知，红色地方所受的应力最大，最大应力为：23.76Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.002206μl（5）网格数68788根据应力云图可知，红色地方所受的应力最大，最大应力为：25.97Mpa 根据应变云图可知，红色地方的应变最大，最大应变为：0.002454μl综合上述5种请况可知随着网格的细分，所得的应变以及应力的结果是收敛的。
计算力学(有限元)matlab编程大作业(空间网架)
逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 输入结构约束总数：11 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：3 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：3 按照编号逐一输入结构约束编号（先输入节点号） ） ：4 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：4 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：3 输入施加在结构上的总的力的个数：3 按照编号逐一输入力作用节点编号（先输入节点号） ） ：2 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：-1 输入力的大小：1000 按照编号逐一输入力作用节点编号（先输入节点号） ） ：2 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：1
有限元分析大作业
基于ＡNSYS软件的有限元分析报告机制1２05班杜星宇U201２10671一、概述本次大作业主要利用ANSYS软件对桌子的应力和应变进行分析,计算出桌子的最大应力和应变.然后与实际情况进行比较，证明分析的正确性,从而为桌子的优化分析提供了充分的理论依据,并且通过对ANSYS软件的实际操作深刻体会有限元分析方法的基本思想,对有限元分析方法的实际应用有一个大致的认识。
二、问题分析已知：桌子几何尺寸如图所示,单位为mm。
假设桌子的四只脚同地面完全固定,桌子上存放物品,物品产生的均匀分布压力作用在桌面,压力大小等于３00Ｐa，其中弹性模量E＝9。
3GPa,泊松比μ=０。
３5,密度ρ=560kg/ｍ3，分析桌子的变形和应力．将桌脚固定在地面，然后在桌面施加均匀分布的压力，可以看作对进行平面应力分析，桌脚类似于梁单元。
由于所分析的结构比较规整且为实体，所以可以将单元类型设为八节点六面体单元。
操作步骤如下：1、定义工作文件名和工作标题（1)定义工作文件名:执行UｔｉｌitｙＭenu/ File/ChangｅＪoｂnａme,在弹出Chaｎｇe Ｊobnａme 对话框修改文件名为Tａblｅ。
选择Ｎew log aｎｄerrｏｒｆiｌeｓ复选框。
(2)定义工作标题：Ｕtility Ｍenu／File／Ｃhaｎｇe Tｉtle，将弹出ChangｅTit ｌｅ对话框修改工作标题名为Tｈｅ aｎａlysiｓ of taｂｌe。
（3)点击:Pｌｏt/Ｒｅplot。
2、设置计算类型（１）点击:Main Menu/Preferencｅs，选择Structurａl，点击OK。
3、定义单元类型和材料属性(1)点击:Ｍain Meｎu／Prｅprｏcｅssｏr／Ｅlement Type/Add/Edit/Deleｔｅ,点击Aｄd，选择Solid〉Brｉｃk ８ｎode 185,点击OK,点击Clｏse。
（2)点击Main ｍenu/preｐｒｏceｓsｏr/Matｅｒiaｌ Prｏps/Maｔerial Models / Struｃtural／ Liｎｅar/ Elａstiｃ／Ｉsotropic，设置EＸ为9.３e9,PRXY为0。
- 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个结点上的应力可由高斯点上的应力外推得到的原理，可能会存在一定的误差。