有限元程序设计讲解版

有限单元法程序设计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）结构分析的基本理论，特别是用有限元法对结构进行分析的原理、方法和步骤。

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

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

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

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

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

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

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

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

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