有限单元法-作业..doc

题目： 如图1-1所示，在不考虑塑性变形的条件下视该梁为一悬臂弹性深梁，其中高度，梁长，厚度，外加均布载荷。相关参数为 弹性模量 弹性深梁的有限元模型 　弹性深梁划分网格示意图 为简化计算，将原有的三维问题的悬臂梁视为二维问题。 将该悬臂梁按照“三结点三角形”等参单元形式划分，坐标采用广义坐标（）形式，考虑材料参数、几何尺寸、边界约束以及计算精度等因素影响，将该梁按照85个结点划分为128个单元。 程序设计 一、流程图 二、程序变量说明 前处理程序中的变量 XY(85,2) 将弹性深梁划分为85个结点，2表示各结点的x,y方向唯一 IJM(128,3) 划分为深梁之后产生的128个单元号，3表示单元的结点 有限元主体程序中的变量 B(3,6) 单元刚度矩阵，共计3个节点，6个 BT(6,3) B的转置矩阵 BT1(6,3) 作为中间变量乘积，用于BT与D的乘积存储 D(3,3) 弹性矩阵3×3 P(170) 表示85个结点的受力情况 PXZ（170）对P的修正，为一170×1的列阵 A(170) 85个结点的x,y方向位移 AE(6) 每个单元的各结点x,y方向位移 CJZ(3,6) 用于BT1与B的乘积存放 SIGMA(3) 初始计算后的128个单元的应力。 SIGMA2(3) 应力修匀后的应力值，其中（3）表示的数值 ZJJZ(128,3) 将初始计算出的128个单元值存储，以便下一步的计算数值 KE(6,6)　单元刚度矩阵 K(170,170)　深梁总体刚度矩阵，共170×170个元素 KXZ(170,170) 修正后的深梁总体刚度矩阵，元素数与K一样 KNZ(170,170) 对K矩阵求逆后的矩阵， 三、总程序 C PROGRAM DIMENSION XY(85,2),IJM(128,3),B(3,6),BT(6,3),BT1(6,3),D(3,3),P(170 *),PXZ(170),A(170),AE(6),CJZ(3,6),SIGMA(3),SIGMA2(3),ZJJZ(128,3) REAL KE(6,6),K(170,170),KXZ(170,170),KNZ(170,170) OPEN(1,FILE=JZJD.DAT,STATUS=UNKNOWN) DO I=1,16 DO J=1,4 MS=(I-1)*8+J MX=(I-1)*8+J+4 IJM(MS,1)=(I-1)*5+J IJM(MS,2)=IJM(MS,1)+6 IJM(MS,3)=IJM(MS,1)+1 IJM(MX,1)=IJM(MS,1) IJM(MX,2)=IJM(MS,1)+5 IJM(MX,3)=IJM(MS,1)+6 END DO END DO DO I=1,128 WRITE(1,10) I,IJM(I,1),IJM(I,2),IJM(I,3) END DO 10 FORMAT(1X,6I3) C PAUSE DO I=1,17 DO J=1,5 N=(I-1)*5+J XY(N,1)=(I-1)*0.25 XY(N,2)=(J-1)*0.25 WRITE(1,*) N,XY(N,1),XY(N,2) END DO END DO C PAUSE CLOSE(1) OPEN(2,FILE=DGJZ.DAT,STATUS=UNKNOWN) OPEN(1,FILE=JZJD.DAT,STATUS=UNKNOWN) OPEN(3,FILE=BJZ.DAT,STATUS=UNKNOWN) E=210*1.E9 U=0.3 DO I=1,3 DO J=1,3 D(1,1)=E/(1-U**2) D(1,2)=U*E/(1-U**2) D(1,3)=0 D(2,1)=U*E/(1-U**2) D(2,2)=E/(1-U**2) D(2,3)=0 D(3,1)=0 D(3,2)=0 D(3,3)=E*(1-U)/2*(1-U**2) C WRITE(*,*) D(I,J) END DO END DO !计算B的数值，遍历128次 A1=1./32 DO N=1,128 B(1,1)=1./(2*A1)*(XY(IJM(N,2