有限元分析程序设计

有限元分析程序设计学校：燕山大学院系：建筑工程与力学学院专业：01级工程力学姓名：张任良指导老师：杜国君完成时间2004年12月31日连续体平面问题的有限元程序分析[题目]：如下图的正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力，，同时在沿对角线y轴上受一对集中压力，载荷为2KN，假设取板厚，泊松比。2kN2kN2kN1kN/m[分析过程]：由于连续平板的对称性，只需要取其在第一象限的四分之一局部参加分析，然后人为作出一些辅助线将平板“分割”成假设干局部，再为每个局部选择分析单元。采用将此模型化分为4个全等的直角三角型单元。利用其对称性，四分之一局部的边界约束，载荷可等效如下图。1kN/m1kN/m[程序原理及实现]：用FORTRAN程序的实现。由节点信息文件NODE.IN和单元信息文件ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。模型根本信息由文件为BASIC.IN生成。该程序的特点如下：问题类型：可用于计算弹性力学平面问题和平面应变问题单元类型：采用常应变三角形单元位移模式：用用线性位移模式载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷材料性质：弹性体由单一的均匀材料组成约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束方程求解：针对半带宽刚度方程的Gauss消元法输入文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT.IN结果文件：输出一般的结果文件DATA.OUT程序的原理如框图：开始开始输入数据〔子程序READ_IN〕BASIC.IN〔根本信息文件〕NODE.IN〔节点信息文件〕ELEMENT.IN(单元信息文件)形成单元刚度矩阵〔子程序FORM_KE〕以半带存储方式形成整体刚度矩阵〔BAND_K〕形成节点载荷向量〔子程序FORM_P〕处理边界条件〔子程序DO_BC〕求解方程获得节点位移〔子程序SOLVE〕计算单元及节点应力〔子程序〕结束输出方件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坐标值LJK_U(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的格式如下表栏目格式说明实际需输入的数据根本模型数据第1行，每两个数之间用“，”号隔开问题类型，单元个数，节点个数，有约束的节点数，有载何的节点数材料性质第2行，每两个数之间用“，”号隔开弹性模量，泊松比，单元厚度节点约束信息在材料性质输入行之后另起行，每两个数之间用“，”号隔开LJK_U(N_BC,3)位移约束的节点编号，该节点x方向约束代码，该节点y方向代码，节点荷载信息在节点约束信息输入行之后另起行，每两个数之间用“，”号隔开P_IJK(N_LOAD,3)载荷作用的节点编号，该节点x主向载荷，该节点y方向载荷，……需读入的节点信息文件NODE.IN的格式如下表栏目格式说明实际需输入的数据节点信息每行为一个节点的信息〔每行三个数，每两个数之间用空格或“，”分开〕ND_ANSYS(N_NIDE)节点号，该节点的x坐标，该节点y方向坐标需读入的单元信息文件ELEMENT.IN的格式如下表栏目格式说明实际需输入的数据单元信息每行为一个单元的信息〔每行有14个整型数，前4个为单元节点编号，对于3节点编号，第4个节点编号与第3个节点编号相同，后10个数无用，可输入“0”，每两个整型数之间用至少一个空格分开〕NE_ANSYS(N_ELE,14)单元的节点号1〔空格〕单元的节点号2〔空格〕单元的节点号3〔空格〕单元的节点号4〔空格〕0〔空格〕0〔空格〕0〔空格〕0〔空格〕0〔空格〕0〔空格〕0〔空格〕0〔空格〕0〔空格〕0输出结果文件DATA.OUT格式如下表栏目实际输出的数据节点位移IRESULT_N(2*I_1)RESULT_N(2*I)节点号x方向位移y方向位移单元应力的三个分量IESTE_ELE(IE,1)STE_ELE(IE,2)STE_ELE(IE,3)单元号x方向应力y方向应力剪切应力节点应力的三个分量ISTS-ND(I,1)STS-ND(I,2)STS-ND(I,3)节点号x方向应力y方向应力剪切应力[算例原始数据和程序分析]：〔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的数据为10.02.020.01.031.01.040.0.51.00.62.00.〔3〕手工准备的单元信息文件ELEMENT.IN的数据为12330000111101245500001111025322000011110335660000111104〔4〕源程序文件chengxu.for为：PROGRAMFEM2D DIMENSIONIJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),&RESULT_N(500),AK(500,100) DIMENSIONSTS_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_LOAD IF(ID.EQ.1)WRITE(8,20) IF(ID.EQ.2)WRITE(8,25)20 FORMAT(/5X,'=========PLANESTRESSPROBLEM========')25 FORMAT(/5X,'=========PLANESTRAINPROBLEM========') CALLREAD_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT,&IJK_ELE,X,Y,IJK_U,P_IJK) CALLBAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,& IJK_ELE,X,Y,PE,PR,PT,AK) CALLFORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,& RESULT_N) CALLDO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) CALLSOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) CALLCAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N,&STS_ELE,STS_ND)ctoputoutadatafile WRITE(9,70)REAL(N_NODE),REAL(N_ELE)70FORMAT(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)71FORMAT(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) ENDc ctogettheoriginaldatainordertomodeltheproblemSUBROUTINEREAD_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,&PT,IJK_ELE,X,Y,IJK_U,P_IJK)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3),& P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14) REALND_ANSYS(N_NODE,3) READ(4,*)PE,PR,PT READ(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) DO10I=1,N_NODE X(I)=ND_ANSYS(I,2) Y(I)=ND_ANSYS(I,3)10CONTINUEDO11I=1,N_ELE DO11J=1,3IJK_ELE(I,J)=NE_ANSYS(I,J)11CONTINUE N_BAND=0 DO20IE=1,N_ELE DO20I=1,3 DO20J=1,3 IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J)) IF(N_BAND.LT.IW)N_BAND=IW20 CONTINUEN_BAND=(N_BAND+1)*2 IF(ID.EQ.1)THEN ELSE PE=PE/(1.0-PR*PR) PR=PR/(1.0-PR) ENDIF RETURN ENDcCtoformthestiffnessmatrixofelementSUBROUTINEFORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3),& AKE(6,6),SS(6,6) CALLCAL_DD(PE,PR,DD) CALLCAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DO10I=1,3 DO10J=1,6 SS(I,J)=0.0 DO10K=1,310SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO20I=1,6 DO20J=1,6 AKE(I,J)=0.0 DO20K=1,320AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PT RETURN ENDcctoformbandedglobalstiffnessmatrix SUBROUTINEBAND_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_NODE DO40I=1,N_DOF DO40J=1,N_BAND40 AK(I,J)=0DO50IE=1,N_ELE CALLFORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DO50I=1,3 DO50II=1,2 IH=2*(I-1)+II IDH=2*(IJK_ELE(IE,I)-1)+II DO50J=1,3 DO50JJ=1,2 IL=2*(J-1)+JJ IZL=2*(IJK_ELE(IE,J)-1)+JJ IDL=IZL-IDH+1 IF(IDL.LE.0)THEN ELSE AK(IDH,IDL)=AK(IDH,IDL)+AKE(IH, ENDIF50 CONTINUE RETURN ENDcctocalculatetheareaofelementSUBROUTINECAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE) DIMENSIONIJK_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.0 RETURN ENDcctocalculatetheelasticmatrixofelement SUBROUTINECAL_DD(PE,PR,DD) DIMENSIONDD(3,3) DO10I=1,3 DO10J=1,310DD(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) RETURN ENDc ctocalculatethestrain-displacementmatrixofelementSUBROUTINECAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DIMENSIONIJK_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) DO10II=1,3 DO10JJ=1,310BB(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) CALLCAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE) DO20I1=1,3 DO20J1=1,620 BB(I1,J1)=BB(I1,J1)/(2.0*AE)RETURN ENDcctoformthegloballoadmatrix SUBROUTINEFORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,& RESULT_N) DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3),& RESULT_N(N_DOF) DO10I=1,N_DOF10 RESULT_N(I)=0.0 DO20I=1,N_LOAD II=P_IJK(I,1) RESULT_N(2*II-1)=P_IJK(I,2)20 RESULT_N(2*II)=P_IJK(I,3)RETURNENDcctodealwithBC(u)(hereonlyforfixeddisplacement)using"1-0"method SUBROUTINEDO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) DIMENSIONRESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100) DO30I=1,N_BC IR=IJK_U(I,1) DO30J=2,3 IF(IJK_U(I,J).EQ.0)THEN ELSE II=2*IR+J-3 AK(II,1)=1.0 RESULT_N(II)=0.0 DO10JJ=2,N_BAND10AK(II,JJ)=0.0DO20JJ=2,II20AK(II-JJ+1,JJ)=0.0ENDIF30CONTINUERETURN ENDcctosolvethebandedFEMequationbyGAUSSelimination SUBROUTINESOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) DIMENSIONRESULT_N(N_DOF),AK(500,100) DO20K=1,N_DOF-1 IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1 IF(N_DOF.LE.K+N_BAND-1)IM=N_DOF DO20I=K+1,IM L=I-K+1 C=AK(K,L)/AK(K,1)IW=N_BAND-L+1 DO10J=1,IW M=J+I-K10AK(I,J)=AK(I,J)-C*AK(K,M)20RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K)RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1) DO40I1=1,N_DOF-1 I=N_DOF-I1 IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1 IF(N_BAND.LE.N_DOF-I-1)JQ=N_BAND DO30J=2,JQ K=J+I-130RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)40RESULT_N(I)=RESULT_N(I)/AK(I,1)WRITE(8,50)50 FORMAT(/12X,'*****RESULTSBYFEM2D*****',//8X,&'--DISPLACEMENTOFNODE--'//5X,'NODENO',8X,'X-DISP',8X,'Y-DISP') DO60I=1,N_NODE60WRITE(8,70)I,RESULT_N(2*I-1),RESULT_N(2*I)70FORMAT(8X,I5,7X,2E15.6)RETURN ENDcccalculatethestresscomponentsofelementandnode SUBROUTINECAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N,&STS_ELE,STS_ND) DIMENSIONIJK_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)DIMENSIONSTS_ELE(500,3),STS_ND(500,3) WRITE(8,10)10FORMAT(//8X,'--STRESSESOFELEMENT--')CALLCAL_DD(PE,PR,DD) DO50IE=1,N_ELE CALLCAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DO20I=1,3 DO20J=1,6 SS(I,J)=0.0 DO20K=1,320SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO30I=1,3 DO30J=1,2 IH=2*(I-1)+J IW=2*(IJK_ELE(IE,I)-1)+J30DISP_E(IH)=RESULT_N(IW)STX=0 STY=0 TXY=0 DO40J=1,6STX=STX+SS(1,J)*DISP_E(J)STY=STY+SS(2,J)*DISP_E(J)40TXY=TXY+SS(3,J)*DISP_E(J)STS_ELE(IE,1)=STXSTS_ELE(IE,2)=STY STS_ELE(IE,3)=TXY50WRITE(8,60)IE,STX,STY,TXY60FORMAT(1X,'ELEMENTNO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)cthefollowingpartistocalculatestresscomponentsofnode WRITE(8,55)55FORMAT(//8X,'--STRESSESOFNODE--')DO90I=1,N_NODE A=0. B=0. C=0. II=0 DO70K=1,N_ELE DO70J=1,3 IF(IJK_ELE(K,J).EQ.I)THEN II=II+1A=A+STS_ELE(K,1)B=B+STS_ELE(K,2) C=C+STS_ELE(K,3) ENDIF70CONTINUESTS_ND(I,1)=A/II STS_ND(I,2)=B/II STS_ND(I,3)=C/II WRITE(8,75)I,STS_ND(I,1),STS_ND(I,2),STS_ND(I,3)75FORMAT(1X,'NODENO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)90CONTINUERETURN ENDcFEM2Dprogrammend[算例结果]：chengxu.for所输出的数据文件DATA.OUT数据内容如下：=========PLANESTRESSPROBLEM========*****RESULTSBYFEM2D*****--DISPLACEMENTOFNODE--NODENOX-DISPY-DISP1.000000E+00-.525275E+012.000000E+00-.225275E+013-.108791E+01-.137363E+014