有限元程序设计讲解版

有限单元法程序设计123456789101112( 1 )( 2 )( 3 )( 4 )( 5 )( 6 )( 7 )( 8 )( 9 )( 10 )( 11 )( 12 )( 13 )( 14 )( 15 )程序一：平面刚架静力分析程序（PF.FOR ）用平面刚架静力分析程序分析所示结构。

图示刚架材料为钢筋混凝土，杆件截面为矩形，柱截面宽0.4m ，高0.4m ；梁截面宽0.4m ，高0.7m 。

材料的弹性模量E=3.55E7。

1． 原始数据的准备与输入20KNm 6 m 6 m 4 m3.5 m3.5 m（4）单元荷载将以上数据以下面的方式输入到程序中：************************************************************* * * * qiuchangting gangjia 2009.6 * * * ************************************************************** 3.55E7 15 12 9 11 2 0.16 2.133E-32 3 0.16 2.133E-33 4 0.16 2.133E-34 5 0.28 11.43E-35 6 0.16 2.133E-36 7 0.16 2.133E-37 8 0.16 2.133E-35 9 0.28 11.43E-39 10 0.16 2.133E-310 11 0.16 2.133E-311 12 0.16 2.133E-33 6 0.28 11.43E-36 10 0.28 11.43E-32 7 0.28 11.43E-37 11 0.28 11.43E-30 00 40 7.50 116 116 7.56 46 012 1112 7.512 412 011 012 013 081 082 083 0121 0122 0123 051 3 20 42 3 20 3.53 3 20 3.513 4 -10 614 4 -10 6程序运行后，输出结果为：************************************************************* * * * qiuchangting gangjia 2009.6 * * * ************************************************************** The Input DataThe General InformationE NM NJ NS NLC3.550E+07 15 12 9 1The Information of Membersmember start end A I1 12 1.600000E-01 2.133000E-032 23 1.600000E-01 2.133000E-033 34 1.600000E-01 2.133000E-034 45 2.800000E-01 1.143000E-025 56 1.600000E-01 2.133000E-036 67 1.600000E-01 2.133000E-037 7 8 1.600000E-01 2.133000E-038 5 9 2.800000E-01 1.143000E-029 9 10 1.600000E-01 2.133000E-0310 10 11 1.600000E-01 2.133000E-0311 11 12 1.600000E-01 2.133000E-0312 3 6 2.800000E-01 1.143000E-0213 6 10 2.800000E-01 1.143000E-0214 2 7 2.800000E-01 1.143000E-0215 7 11 2.800000E-01 1.143000E-02The Joint Coordinatesjoint X Y1 .000000 .0000002 .000000 4.0000003 .000000 7.5000004 .000000 11.0000005 6.000000 11.0000006 6.000000 7.5000007 6.000000 4.0000008 6.000000 .0000009 12.000000 11.00000010 12.000000 7.50000011 12.000000 4.00000012 12.000000 .000000The Information of SupportsIS VS11 .00000012 .00000013 .00000081 .00000082 .00000083 .000000121 .000000122 .000000123 .000000( NA= 306 )( NW= 1038 )Loading Case 1The Loadings at JointsNLJ= 0The Loadings at MembersNLM= 5ILM ITL PV DST1 3 20.0000 4.0000002 3 20.0000 3.5000003 3 20.0000 3.50000013 4 -10.0000 6.00000014 4 -10.0000 6.000000The Results of CalculationThe Joint Displacementsjoint u v phi1 5.517799E-21 3.746750E-21 -1.212291E-202 5.035124E-03 2.638557E-05 -5.743715E-043 7.666356E-03 4.137788E-05 -2.016188E-044 8.569997E-03 4.229713E-05 -1.643039E-055 8.555549E-03 -5.769378E-05 -3.895922E-056 7.633712E-03 -6.086508E-05 -1.701734E-047 5.006871E-03 -4.326034E-05 -8.670744E-058 6.862436E-21 -6.142968E-21 -1.388901E-209 8.548212E-03 -1.060822E-04 -7.533107E-0510 7.621811E-03 -1.019917E-04 -1.262908E-0411 4.992188E-03 -6.763227E-05 -5.169941E-0412 5.619765E-21 -9.603783E-21 -1.221822E-20The Terminal Forcesmember N(st) Q(st) M(st) N(en) Q(en) M(en)1 -37.468 95.178 147.896 37.468 -15.178 72.8162 -24.330 61.984 59.575 24.330 8.016 34.8703 -1.492 46.064 35.772 1.492 23.936 2.9524 23.936 -1.492 -2.952 -23.936 1.492 -5.9995 -5.147 11.780 23.454 5.147 -11.780 17.7776 28.570 46.144 78.946 -28.570 -46.144 82.5587 61.430 68.624 135.607 -61.430 -68.624 138.8908 12.156 -6.638 -17.455 -12.156 6.638 -22.3759 6.638 12.156 22.375 -6.638 -12.156 20.17010 55.760 31.872 64.229 -55.760 -31.872 47.32311 96.038 56.198 102.608 -96.038 -56.198 122.18212 54.080 -22.839 -70.642 -54.080 22.839 -66.38913 19.716 10.878 -30.334 -19.716 49.122 -84.39814 46.806 -13.137 -132.391 -46.806 73.137 -126.43215 24.326 -40.277 -91.733 -24.326 40.277 -149.931( NA= 306 )( NW= 1058 )源程序：C PF.FOR (A program for analysis of plane frame)C Version 4.3 1994C Main program reads the control data & organizes the wholeC calculation by calling subroutines.DIMENSION W(20000)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)C Read data of members, joints, supports & print themDIMENSION 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)C 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) C Determine location of diagonal elements of global stiffnessC 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)C Calculate length, cosine & sine of memberDIMENSION 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)C Calculate element stiffness matrix along local axesDIMENSION 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)C Calculate element stiffness matrix along global axesDIMENSION 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)C 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)C Introduce support conditions into global stiffness matrix ADIMENSION 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)C Solve equations (1) - decomposition of matrix ADIMENSION 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) C Solve equations (2) - forward & back substitutionDIMENSION 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)C 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)C Read data of loading at joint & form joint load vector BDIMENSION 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 FORMA T(/2X,'ILJ',11X,'PX',14X,'PY',15X,'PM')WRITE (*,2)(ILJ(M),PX(M),PY(M),PM(M),M=1,NLJ)2 FORMA T(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)C Initialize terminal forces of membersDIMENSION F(6,NM)WRITE (*,1) NLM1 FORMA T(/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)C Read data of loading at member & calculate fixed-end forces,C 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)*SEND IFIF (ITL(M).EQ.2.OR.ITL(M).EQ.4) THENP1=PV(M)*SP2=PV(M)*CEND IFIF (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)END IFIF (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-F4END IFF(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(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)C Introduce support conditions into joint load vector BDIMENSION 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)C 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 FORMA T(1X,I4,1P3E16.6)RETURNENDSUBROUTINE COTF (E,NM,NJ,N,IST,IEN,AR,RI,X,Y,LV,B,F)C Calculate & print terminal forces of membersDIMENSION 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 FORMA T(/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))U5=-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程序二：平面三节点有限元程序如下图所示的悬臂梁，受均布荷载q=1N/mm2 作用。

结构有限元分析程序设计绪论§0.1 开设“有限元程序设计”课程的意义和目的§0.2 课程特点§0.3 课程安排§0.4 课程要求§0.5 基本方法复习$0.1 意义和目的1.有限元数值分析技术本身要求工程设计研究人员掌握1). 有限元数值分析技术的完善标志着现代计算力学的真正成熟和实用化,已在各种力学中得到了广泛的应用。

比如：，已杨为工程结构分析中最得以收敛的技术手段，现代功用大致有：a). 现代结构论证。

对结构设计从内力，位移等方面进行优劣评定，从而进行结构优化设计。

b)可取代部份实验，局部实验+有限元分析，是现代工程设计研究方法的一大特点。

c)结构的各种功能分析（疲劳断裂，可靠性分析等）都以有限元分析工具作为核心的计算工具。

2). 有限元数值分析本身包括着理论+技术实现（本身功用所绝定的）有限元数值分析本身包括着泛函理论+分片插值函数+程序设计2. 有限元分析的技术实现（近十佘年的事）更依赖于计算机程序设计有限元分析的技术取得的巨大的成就，从某种意义上说，得益于计算机硬件技术的发展和程序设计技术的发展，这两者的依赖性在当代表现得更加突出。

（如可视化技术）3.从学习的角度，不仅要学习理论，而且要从程序设计设计角度对这些理论的技术实现有一个深入的了解，应当致力于掌握这些技术实现能力，从而开发它，发展它。

（理论本身还有待于进一步完美相应的程序设计必须去开发）4.程序设计不仅是实现有限元数值分析的工具和桥梁，而且在以下诸方面也有意义：1). 精通基本概念，深化理论认识；2). 锻炼实际工程分析，实际动手的能力；3). 获得以后工作中必备的工具。

（作业+老师给元素库）目的：通过讲述有限元程序设计的技术与技巧，便能达到自编自读的能力。

§0.2 课程特点总描述：理论+算法+数据结构（程序设计的意义）理论：有限元算法，构造，步骤，解的等外性，收敛性，稳定性，误差分析算法；指求解过程的技术方法，含两方面的含义；a. 有限元数值分析算法，b, 与数据结构有关的算法（总刚稀疏存贮，提取，节点优化编号等）数据结构：指各向量矩阵存贮管理与实现，辅助管理结构（指针，数据记录等）具体特点：理论性强：能量泛函理论+有限元构造算法+数据结构构造算法内容繁杂：理论方法+技术方法+技术技巧技巧性强：排序，管理结构（指针生成，整型运算等）§0.3 课程安排①. 单元刚度矩阵及元素设计（单元刚阵算法，杆梁平面分析，板弯非协调元等）②. 总刚的形式及程序设计（单刚提前准备，技术复杂）③. l边界条件及程序设计（等效荷载计算，位移边界条件置入，多工况的对称性）④. 总刚线性方程组求解（LDL T分解，分块算法，子结构算法，波前法）⑤．单元应力计算+应力处理与改善。

第9章有限元法程序设计9．1 引言在用有限元法进行结构分析时，将会遇到大量的数值计算，因而在实用上是一定要借助于计算机和有限元程序，才能完成这些复杂而繁重的数值计算工作。

事实上，有限元程序的设计是有限元研究的一个很重要的部分。

它是理论和方法的载体，是理论用于实际必不可少的桥梁，是有限元学术研究与实际应用水平的代表。

有好的、高深的理论和算法并不等于有好的程序，还必须有实际的程序开发经验的多年积累、丰富的计算机知识、大量的资金和人力的投入，多年的开发修正与改进才能编制出好的程序来。

一些著名的有限元程序开发的发展历史也体现出了这一规律。

设计一个用于结构分析的有限元法程序，要求设计者至少应该掌握下列知识：（１）掌握一种程序开发工具，如VC（Visual C++），CB（C++Buildel），Delphi，VB（Visual Basic）或VF（Visual Fortran）等。

在本书中所有程序均用VC写出。

（2）数值方法，如线性和非线性代数方程的求解，矩阵特征值的求解以及数值积分等。

（3）结构分析的基本理论，特别是用有限元法对结构进行分析的原理、方法和步骤。

由于一般的软件工程师不懂结构分析原理，因此，结构分析程序的开发任务主要应由结构工程师来承担。

掌握结构分析程序设计方法，是以计算机辅助设计为主要标志的现代工程设计方法对结构工程师的要求。

作为结构工程师，应该具有对结构分析程序的使用、阅读、修改和编制的基础知识和技术素质。

有限元程序的总体组成可分为三个部分：前处理部分，有限元分析本体部分以及后处理部分。

有限元分析本体部分是有限元分析程序的核心。

它根据离散模型的数据文件进行有限元分析，有限元分析的原理和采用的数值方法集中于此。

因此，这一部分程序是有限元分析是否准确可靠的关键部分。

有限元分析所使用的离散模型的数据文件主要包括：模型的节点数、节点坐标与节点编码，单元数据与单元编码；材料和载荷信息等。

实际工程问题的离散模型数据文件十分庞大。

有限单元法及程序设计绪论1.力学分析方法:解析法,数值法有限元法——实际结构形状和所受载荷比较复杂,大多用解析法很困难,因而数值法得到不断发展,随着电子计算机的进步,而发展起来的一种新兴的数值分析方法.2基本步骤：（1）结构离散化：将结构从集合上用线或面划分为有限个单元。

（2）单元分析：导出单元的节点位移和结点力之间的关系（单元刚度矩阵）。

（3）整体分析：将各单元组成的结构整体进行分析，导出征个结构点位移与结点力之间的关系。

3程序设计的步骤：（1）提出问题，拟定解决方案（2）构造数学模型（3）画出程序流程图（4）编写程序（5）编译调试程序（6）试算验证程序4.根据国家标准（GB-1526-89）规定的程序流程图标准化符号及规定：a)图表示程序流程图的起点和终点；b)图表示数据信息的输入和输出；c)图表示数据进行系列运算之前要完成的数据预置；d)图表示判断条件；e)图表示各种处理功能，如数学运算方式等；f)图表示流程的路径和指向。

第一篇杆件结构的有限单元法及程序设计第一章平面杆件单元的有限单元法第一节有限单元法的基本概念1．基本思路：先分后合（先单元分析，再整体分析）2．基本概念：整体号：节点端点号按自然数1，2，3，……（在整体坐标系xOy下）局部号：每一个单元始末用i，j标记（在单元的局部坐标xyz系下，方向与整体坐标系一致）。

⎡k k⎤ii ij 3．F e=k eδe其中：ke=⎥⎢k kji jj⎣⎦单元刚度矩阵，各元素为刚度系数⎡θ⎤iδe=⎥⎢θ⎣j⎦单元杆段位移列阵⎡M⎤iF e=⎥⎢Mj⎣⎦单元杆端力列阵K∆=P（1-7）K=⎡k112131⎢k⎢⎢k⎣k12k22k32k⎤13⎥k23⎥k⎥33k⎥⎦整体刚度矩阵∆=[]θ1θθT23位移列阵P=[]M1M M节点载荷列阵233.有限元位移法分析连续梁需要考虑的问题(1)刚度集成法:①将(1-3)K扩阶,扩大的元素为0,得到单元贡献矩阵单元①:K①=1⎡kii⎢k⎢ji⎢02kijkjj310⎤1⎡0⎥⎢020⎥⎢单元②:K=②2kiikji30⎤⎥kij⎥k⎥123⎣0⎥3⎢0⎦⎣jj ⎦②将单元贡献矩阵想叠加,形成整体刚度矩阵123⎡k k0⎤11ii ij⎢⎥K=K①+K②=k k k k 2⎢⎥1122ji jj+ii ji⎢0k2⎥32k⎣⎦ji jj(2)两端支承条件的引入先不考虑约束条件,得到整体刚度矩阵后,将其主对角线元素k ii改为1,第i行,第j列其余元素改为0,对应的载荷元素也改为0.(3)非结点荷载的处理利用等效结点荷载进行分析:1各结点(包括两端结点)加约束,阻止结点转动,其约束力矩分别为交于该结点的各相关单元的固端力矩之和,顺时针为正.2去掉附加约束(相当在各结点施加外力荷载P3将两部分杆端弯矩叠加起来.e,其大小与约束力矩相同,方向相反)第二节局部坐标系中的单元刚度矩阵1.一般单元设单元○e的弹性模量、截面惯性矩、截面积分别为E、I、A,杆长为l。

有限元法基础讲稿-第3讲新1.12.23.3将第步生成的文件复制到目录下如此目录中已有该文件则覆盖它，问题描述如图所示为一个悬臂梁在力作用下求该梁点的挠度，如图所示输入数据然后输入及，弹性体在载荷作用下体内任意一点的应力状态可由个应力分量。

有限元法基础讲稿-第3讲新2017-08-16 09:26:07 | #1楼ANSYS 基本操作INTRODUCTION TO ANSYS 5.7 - Part 1ANSYS界面与操作无论版本怎样变化，始终以原貌为主，仅作少量的改进，具有较强的继承性，形成了自己固有的风格，具有操作直观易行的特点。

论述以符号“>”表示进入下一级菜单或选择项。

主题:A. ANSYS安装B. ANSYS启动、用户界面及退出C. ANSYS操作方式D. ANSYS典型分析过程A. ANSYS安装INTRODUCTION TO ANSYS 5.7 - Part 1找出你的计算机名及MAC地址。

ANSYS安装光盘放入后会自动运行，然后点击最后一行“Display the license server hostid”，会弹出一个窗口，其中第一行 HOSTNAME是计算机名，第二行FLEXID 为MAC地址！编辑“ansys.dat”文件。

把光盘内的MAGNiTUDE(或有“ansys.dat”文件的目录)目录复制一份到硬盘上，然后用记事本编辑“ansys.dat”文件。

第一行：其中“host”用计算机名代替，“00000000000”用MAC地址（12位）代替。

然后保存。

生成“license.dat”。

运行“keygen.bat”，即可立即生成一个“license.dat”文件。

安装ANSYS10.0。

运行“setup.exe”，全部按默认点击即可。

当出现“Is this a license server machine?”时点“是”。

出现“ANSYS FLEXlm license file”对话框时，寻Browse for the location ofan existing license file”选项，然后指向生成的“license.dat”文件和“ansys.dat”文件。

有限元教材-第十章有限元程序设计第十章有限元程序设计有限元方法作为一门系统的技术，仅学会了它的基本理论是远远不够的，只有形成完整的计算程序，问题才最终得到了解决。

完成这样的有限元程序设计是一项工作量很大的工程。

本章就是要结合简单的有限元教学程序FEMED，简要介绍有限元程序设计技术。

FEMED 是专为有限元程序设计教学编制的程序，它不包含复杂的前后处理功能，可进行平面问题及平面桁架的线弹性静力分析，在程序结构上与大型程序类似，具有计算单元的任意扩充功能，在方程的组集和求解上也采用了较为流行的变带宽存储方式。

有限元程序大致可分为两类，第一类是专用程序，主要用于研究或教学，一般这类程序规模较小，前后处理功能较弱。

用于研究的程序能够解一些特殊的问题，满足研究工作的需要。

而教学程序则是为了学生了解有限元的主要结构和设计方法设计的，程序比较简单，FEMED就属于这类程序。

第二类是大型通用程序，是大型结构分析的得力工具，目前国际上流行的大约有2000多种。

常用的有NASTRAN、MARC、ANSYS、ADINA和ABAQUS等。

这类程序一般前后处理功能比较强，有友好的界面，能进行大型计算，但往往无法完成具有特殊要求的计算。

通过本章的学习，使读者初步掌握有限元编程的基本方法，具有开发特殊功能的专用程序或为通用程序开发具有特殊功能的计算模块的能力。

§10.1有限元程序的基本结构有限元程序一般包括三项基本内容：前处理、结构分析和后处理。

早期有限元分析软件的研究重点在于推导新的高效率求解方法和高精度的单元，随着数值分析方法的逐步完善，尤其是计算机内存和运算速度的飞速发展，整个计算系统用于求解运算的时间越来越少，加之求解问题的日益大型化和复杂化，使得数据准备和运算结果的表现问题日益突出。

因此目前几乎所有的商业化有限元程序系统都有功能很强的前后处理模块，这直接关系到分析软件的可推广性。

它是商用有限元软件不可或缺的部分，但它不是有限元的中心部分，在本书中不作详细介绍。