有限单元法程序设计

有限单元法程序设计
12345678910
11
12( 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 m
3.5 m
3.5 m
（4）单元荷载
将以上数据以下面的方式输入到程序中：
************************************************************* * * * qiuchangting gangjia 2009.6 * * * ************************************************************** 3.55E7 15 12 9 1
1 2 0.16 2.133E-3
2 3 0.16 2.133E-3
3 4 0.16 2.133E-3
4 5 0.28 11.43E-3
5 6 0.16 2.133E-3
6 7 0.16 2.133E-3
7 8 0.16 2.133E-3
5 9 0.28 11.43E-3
9 10 0.16 2.133E-3
10 11 0.16 2.133E-3
11 12 0.16 2.133E-3
3 6 0.28 11.43E-3
6 10 0.28 11.43E-3
2 7 0.28 11.43E-3
7 11 0.28 11.43E-3
0 0
0 4
0 7.5
0 11
6 11
6 7.5
6 4
6 0
12 11
12 7.5
12 4
12 0
11 0
12 0
13 0
81 0
82 0
83 0
121 0
122 0
123 0
5
1 3 20 4
2 3 20 3.5
3 3 20 3.5
13 4 -10 6
14 4 -10 6
程序运行后，输出结果为：
************************************************************* * * * qiuchangting gangjia 2009.6 * * * ************************************************************** The Input Data
The General Information
E NM NJ NS NLC
3.550E+07 15 12 9 1
The Information of Members
member start end A I
1 1
2 1.600000E-01 2.133000E-03
2 2
3 1.600000E-01 2.133000E-03
3 3
4 1.600000E-01 2.133000E-03
4 4
5 2.800000E-01 1.143000E-02
5 5
6 1.600000E-01 2.133000E-03
6 6
7 1.600000E-01 2.133000E-03
7 7 8 1.600000E-01 2.133000E-03
8 5 9 2.800000E-01 1.143000E-02
9 9 10 1.600000E-01 2.133000E-03
10 10 11 1.600000E-01 2.133000E-03
11 11 12 1.600000E-01 2.133000E-03
12 3 6 2.800000E-01 1.143000E-02
13 6 10 2.800000E-01 1.143000E-02
14 2 7 2.800000E-01 1.143000E-02
15 7 11 2.800000E-01 1.143000E-02
The Joint Coordinates
joint X Y
1 .000000 .000000
2 .000000 4.000000
3 .000000 7.500000
4 .000000 11.000000
5 6.000000 11.000000
6 6.000000 7.500000
7 6.000000 4.000000
8 6.000000 .000000
9 12.000000 11.000000
10 12.000000 7.500000
11 12.000000 4.000000
12 12.000000 .000000
The Information of Supports
IS VS
11 .000000
12 .000000
13 .000000
81 .000000
82 .000000
83 .000000
121 .000000
122 .000000
123 .000000
( NA= 306 )
( NW= 1038 )
Loading Case 1
The Loadings at Joints
NLJ= 0
The Loadings at Members
NLM= 5
ILM ITL PV DST
1 3 20.0000 4.000000
2 3 20.0000 3.500000
3 3 20.0000 3.500000
13 4 -10.0000 6.000000
14 4 -10.0000 6.000000
The Results of Calculation
The Joint Displacements
joint u v phi
1 5.517799E-21 3.746750E-21 -1.212291E-20
2 5.035124E-0
3 2.638557E-05 -5.743715E-04
3 7.666356E-03 4.137788E-05 -2.016188E-04
4 8.569997E-03 4.229713E-0
5 -1.643039E-05
5 8.555549E-03 -5.769378E-05 -3.895922E-05
6 7.633712E-03 -6.086508E-05 -1.701734E-04
7 5.006871E-03 -4.326034E-05 -8.670744E-05
8 6.862436E-21 -6.142968E-21 -1.388901E-20
9 8.548212E-03 -1.060822E-04 -7.533107E-05
10 7.621811E-03 -1.019917E-04 -1.262908E-04
11 4.992188E-03 -6.763227E-05 -5.169941E-04
12 5.619765E-21 -9.603783E-21 -1.221822E-20
The Terminal Forces
member N(st) Q(st) M(st) N(en) Q(en) M(en)
1 -37.468 95.178 147.896 37.468 -15.178 72.816
2 -24.330 61.984 59.575 24.330 8.016 34.870
3 -1.492 46.06
4 35.772 1.492 23.936 2.952
4 23.936 -1.492 -2.952 -23.936 1.492 -5.999
5 -5.147 11.780 23.454 5.147 -11.780 17.777
6 28.570 46.144 78.946 -28.570 -46.144 82.558
7 61.430 68.624 135.607 -61.430 -68.624 138.890
8 12.156 -6.638 -17.455 -12.156 6.638 -22.375
9 6.638 12.156 22.375 -6.638 -12.156 20.170
10 55.760 31.872 64.229 -55.760 -31.872 47.323
11 96.038 56.198 102.608 -96.038 -56.198 122.182
12 54.080 -22.839 -70.642 -54.080 22.839 -66.389
13 19.716 10.878 -30.334 -19.716 49.122 -84.398
14 46.806 -13.137 -132.391 -46.806 73.137 -126.432
15 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 1994
C Main program reads the control data & organizes the whole
C calculation by calling subroutines.
DIMENSION W(20000)
CHARACTER IDFN*20,TITLE(5)*72
READ (*,'(A12)')IDFN
OPEN (3,FILE=IDFN,STATUS='OLD')
READ (3,'(A72)')(TITLE(M),M=1,5)
READ (3,*)E,NM,NJ,NS,NLC
L1=1
L2=L1+NM
L3=L2+NM
L4=L3+NM
L11=L4+NM
L12=L11+NJ
L21=L12+NJ
L22=L21+NS
L31=L22+NS
L41=L31+6*NM
CALL 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+N
L52=L51+36
L53=L52+NA*2
L54=L53
L61=L54+N*2
NW=L61+6*NM-1
WRITE (*,1)NA,NW
1 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,NLC
READ (3,*)NLJ
L62=L61+NLJ
L63=L62+NLJ
L64=L63+NLJ
L71=L61
L81=L71+6*NM
CALL 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,*)NLM
L82=L81+NLM
L83=L82+NLM
L84=L83+NLM
CALL 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-1
WRITE (*,1)NA,NW
100 CONTINUE
WRITE (*,'(/)')
END
SUBROUTINE IOMJS (TITLE,E,NM,NJ,NS,NLC,IST,IEN,
& AR,RI,X,Y,IS,VS)
C 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)*72
WRITE (*,'(/)')
WRITE (*,1)(TITLE(M),M=1,5)
1 FORMAT(1X,A72)
WRITE (*,2)E,NM,NJ,NS,NLC
2 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)
RETURN
END
SUBROUTINE LCVCT (NM,IST,IEN,LV,NJ,N)
C Determine location vector of element
DIMENSION IST(NM),IEN(NM),LV(6,NM)
DO 100 M=1,NM
I=IST(M)*3
J=IEN(M)*3
LV(1,M)=I-2
LV(2,M)=I-1
LV(3,M)=I
LV(4,M)=J-2
LV(5,M)=J-1
LV(6,M)=J
100 CONTINUE
N=NJ*3
RETURN
END
SUBROUTINE LCDIA (NM,N,LV,MIN,IBDW,LD,MAXBDW,NA) C Determine location of diagonal elements of global stiffness
C matrix A
DIMENSION LV(6,NM),MIN(N),IBDW(N),LD(N)
DO 100 I=1,N
MIN(I)=I
100 CONTINUE
DO 400 M=1,NM
MINLV=LV(1,M)
DO 200 I=2,6
IF (LV(I,M).LT.MINLV) MINLV=LV(I,M)
200 CONTINUE
DO 300 I=1,6
IF (MINLV.LT.MIN(LV(I,M))) MIN(LV(I,M))=MINLV
300 CONTINUE
400 CONTINUE
MAXBDW=0
DO 500 I=1,N
IBDW(I)=I-MIN(I)+1
IF (IBDW(I).GT.MAXBDW) MAXBDW=IBDW(I)
500 CONTINUE
LD(1)=IBDW(1)
DO 600 I=2,N
LD(I)=LD(I-1)+IBDW(I)
600 CONTINUE
NA=LD(N)
RETURN
END
SUBROUTINE RLCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S)
C 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/RL
S=Y1/RL
RETURN
END
SUBROUTINE KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,
& X,Y,C,S,E1,E2,E3,E4)
C 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)/RL
E2=12.0*E*RI(M)/(RL*RL*RL)
E3=0.5*E2*RL
E4=0.6666667*E3*RL
RETURN
END
SUBROUTINE KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)
C 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*S
A2=(E1-E2)*C*S
A3=E1*S*S+E2*C*C
A4=E3*S
A5=E3*C
A6=E4
AE(1,1)=A1
AE(2,1)=A2
AE(2,2)=A3
AE(3,1)=-A4
AE(3,2)=A5
AE(3,3)=A6
AE(4,1)=-A1
AE(4,2)=-A2
AE(4,3)=A4
AE(4,4)=A1
AE(5,1)=-A2
AE(5,2)=-A3
AE(5,3)=-A5
AE(5,4)=A2
AE(5,5)=A3
AE(6,1)=-A4
AE(6,2)=A5
AE(6,3)=0.5*A6
AE(6,4)=A4
AE(6,5)=-A5
AE(6,6)=A6
RETURN
END
SUBROUTINE FORMA (E,NM,NJ,N,NA,IST,IEN,AR,RI,
& X,Y,LV,AE,LD,A)
C Form global stiffness matrix A
DIMENSION 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,NM
CALL KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE)
DO 200 I=1,6
DO 100 J=1,I
IF (LV(I,M).GE.LV(J,M)) THEN
A(LD(LV(I,M))-LV(I,M)+LV(J,M))
& =A(LD(LV(I,M))-LV(I,M)+LV(J,M))+AE(I,J) ELSE
A(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 IF
100 CONTINUE
200 CONTINUE
300 CONTINUE
RETURN
END
SUBROUTINE AS (NS,N,NA,IS,LD,A)
C Introduce support conditions into global stiffness matrix A
DIMENSION IS(NS),LD(N)
DOUBLE PRECISION A(NA)
DO 100 M=1,NS
I=3*(IS(M)/10)-3+MOD(IS(M),10)
A(LD(I))=1D22
100 CONTINUE
RETURN
END
SUBROUTINE LDLT (N,NA,LD,A,T)
C Solve equations (1) - decomposition of matrix A
DIMENSION LD(N)
DOUBLE PRECISION A(NA),T(N),SUM
DO 300 I=2,N
LDI=LD(I)
I1=I-LDI+LD(I-1)+1
DO 200 J=I1,I-1
LDJ=LD(J)
J1=J-LDJ+LD(J-1)+1
IF(I1.GT.J1) J1=I1
SUM=0.0D0
DO 100 K=J1,J-1
SUM=SUM+T(K)*A(LDJ-J+K)
100 CONTINUE
T(J)=A(LDI-I+J)-SUM
A(LDI-I+J)=T(J)/A(LDJ)
A(LDI)=A(LDI)-T(J)*A(LDI-I+J)
200 CONTINUE
300 CONTINUE
RETURN
END
SUBROUTINE SLVEQ (N,NA,MAXBDW,LD,A,B) C Solve equations (2) - forward & back substitution
DIMENSION LD(N)
DOUBLE PRECISION A(NA),B(N)
DO 200 I=2,N
LDI=LD(I)
I1=I-LDI+LD(I-1)+1
DO 100 J=I1,I-1
B(I)=B(I)-A(LDI-I+J)*B(J)
100 CONTINUE
200 CONTINUE
DO 300 I=1,N
B(I)=B(I)/A(LD(I))
300 CONTINUE
DO 500 I=N-1,1,-1
IMIN=I+MAXBDW
IF(IMIN.GT.N) IMIN=N
DO 400 J=I+1,IMIN
LDJ=LD(J)
J1=J-LDJ+LD(J-1)+1
IF(I.GE.J1) B(I)=B(I)-A(LDJ-J+I)*B(J)
400 CONTINUE
500 CONTINUE
RETURN
END
SUBROUTINE B0 (LC,N,NLJ,B)
C Initialize joint load vector B
DOUBLE PRECISION B(N)
WRITE (*,1)LC,NLJ
1 FORMAT(/15X,'Loading Case',I3
& //10X,'The Loadings at Joints'
& //17X,'NLJ=',I4)
DO 100 I=1,N
B(I)=0.0D0
100 CONTINUE
RETURN
END
SUBROUTINE IOLJB (N,NLJ,ILJ,PX,PY,PM,B)
C 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 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,NLJ
I=ILJ(M)*3
B(I-2)=PX(M)
B(I-1)=PY(M)
B(I)=PM(M)
100 CONTINUE
RETURN
END
SUBROUTINE F0 (NLM,NM,F)
C Initialize terminal forces of members
DIMENSION F(6,NM)
WRITE (*,1) NLM
1 FORMA T(/10X,'The Loadings at Members'
& //17X,'NLM=',I4)
DO 200 J=1,NM
DO 100 I=1,6
F(I,J)=0.0
100 CONTINUE
200 CONTINUE
RETURN
END
SUBROUTINE 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 B
DIMENSION 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,NLM
L=ILM(M)
CALL RLCS (L,NM,NJ,IST,IEN,X,Y,RL,C,S)
D1=DST(M)
D2=RL-D1
IF (ITL(M).EQ.1.OR.ITL(M).EQ.3) THEN
P1=PV(M)*C
P2=-PV(M)*S
END IF
IF (ITL(M).EQ.2.OR.ITL(M).EQ.4) THEN
P1=PV(M)*S
P2=PV(M)*C
END IF
IF (ITL(M).EQ.1.OR.ITL(M).EQ.2) THEN
F1=-P1*D2/RL
F4=-P1-F1
F2=-P2*D2*D2*(RL+2.0*D1)/(RL*RL*RL)
F5=-P2-F2
F3=-P2*D1*D2*D2/(RL*RL)
F6=P2*D1*D1*D2/(RL*RL)
END IF
IF (ITL(M).EQ.3.OR.ITL(M).EQ.4) THEN
G=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-F5
F4=-P1*D1*D1/(2.0*RL)
F1=-P1*D1-F4
END IF
F(1,L)=F(1,L)+F1
F(2,L)=F(2,L)+F2
F(3,L)=F(3,L)+F3
F(4,L)=F(4,L)+F4
F(5,L)=F(5,L)+F5
F(6,L)=F(6,L)+F6
B(LV(1,L))=B(LV(1,L))-F1*C+F2*S
B(LV(2,L))=B(LV(2,L))-F1*S-F2*C
B(LV(3,L))=B(LV(3,L))-F3
B(LV(5,L))=B(LV(5,L))-F4*S-F5*C
B(LV(6,L))=B(LV(6,L))-F6
100 CONTINUE
RETURN
END
SUBROUTINE BS (NS,N,IS,VS,B)
C Introduce support conditions into joint load vector B
DIMENSION IS(NS),VS(NS)
DOUBLE PRECISION B(N)
DO 100 M=1,NS
I=3*(IS(M)/10)-3+MOD(IS(M),10)
B(I)=VS(M)*1D22
100 CONTINUE
RETURN
END
SUBROUTINE OJD (NJ,N,B)
C Print joint displacements
DOUBLE 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)
RETURN
END
SUBROUTINE COTF (E,NM,NJ,N,IST,IEN,AR,RI,X,Y,LV,B,F)
C 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,U6
WRITE (*,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,NM
CALL 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))*S
U2=-B(LV(1,M))*S+B(LV(2,M))*C
U3=B(LV(3,M))
U5=-B(LV(4,M))*S+B(LV(5,M))*C
U6=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 CONTINUE
RETURN
END
程序二：平面三节点有限元程序
如下图所示的悬臂梁，受均布荷载q=1N/mm2 作用。

E=2.1×105 N/mm2 ， ，厚度h=10mm，现用有限元法分析其位移及应力。

=0.3
梁可视为平面应力状态，按图1尺寸划分为均匀的三角形网格，共80个单元，55个节点，坐标轴及单元与节点的编号如图。

将均不荷载分配到相应的节点上，把有约束的节点51、52、53、54、55视作固定铰链，建立如图1所示的离散化计算模型。

计算源程序如下：
PROGRAM MAIN
DIMENSION SK(300,30),EK(12,12),Q(300),MC(55),XY(3,100),XYE(3,4),&
QE(12),NX(4,100)
OPEN(7,FILE='INP.DAT')
REWIND 7
READ(7,10)NF,NE,NN,MB,ND,LE,LS
READ(7,12)E,UM,T
10 FORMAT(7I5)
12 FORMAT(3F15.2)
WRITE(*,600)NF,NE,NN,MB,ND,LE,LS,E,UM,T
ME=NE*NF
MS=NN*NF
CALL INPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)
WRITE(*,102)((XY(I,J),I=1,LS),J=1,NN)
102 FORMAT(10X,'XY'/,(2X,6F12.3))
WRITE(*,101)(Q(I),I=1,MS)
101 FORMAT(10X,'Q'/,(2X,6F12.3))
WRITE(*,500)((NX(I,J),I=1,NE),J=1,LE)
WRITE(*,400)(MC(I),I=1,ND)
500 FORMAT(10X,'NX'/,(2X,12I6))
600 FORMAT(10X,'NF NE NN MB ND LE LS E UM T'/7(2X,I4),3(2X,F8.4)) 400 FORMAT(10X,'MC'/,(2X,10I6))
CALL STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,& UM,T)
CALL SOLVE(SK,Q,MS,MB)
OPEN(9,FILE='OUT.DAT')
REWIND 9
WRITE(9,200)
WRITE(9,250)(Q(I),I=1,MS)
200 FORMAT(5X,'DISPLACEMENT')
250 FORMAT(2X,6E14.5)
CALL STRES(Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T)
STOP 1000
END
SUBROUTINE INPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)
DIMENSION XY(LS,NN),Q(MS),NX(NE,LE),MC(ND)
READ(7,*)XY
READ(7,*)Q
READ(7,*)NX
READ(7,*)MC
CLOSE(7)
10 FORMAT(6F11.2)
20 FORMAT(12I5)
RETURN
END
SUBROUTINE
STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,UM,T) DIMENSION
SK(MS,MB),EK(ME,ME),Q(MS),NX(NE,LE),MC(ND),XY(LS,NN),XYE(LS,NE) DO 35 I=1,MS
DO 35 J=1,MB
35 SK(I,J)=0
DO 200 L=1,LE
DO 40 J=1,NE
LJ=NX(J,L)
DO 40 I=1,LS
40 XYE(I,J)=XY(I,LJ)
DO 50 I=1,ME
DO 50 J=1,ME
50 EK(I,J)=0.0
CALL STIFE(EK,XYE,ME,NE,NF,LS,E,UM,T) IF(L.EQ.1)WRITE(*,70)EK
70 FORMAT(10X,'EK'/,(6E14.5))
DO 200 I=1,NE
DO 200 II=1,NF
M=NF*(I-1)+II
M1=NF*(NX(I,L)-1)+II
DO 200 J=1,NE
DO 200 JJ=1,NF
N=NF*(J-1)+JJ
N1=NF*(NX(J,L)-1)+JJ
MN=N1-M1+1
IF(MN)200,200,150
150 SK(M1,MN)=SK(M1,MN)+EK(M,N)
200 CONTINUE
DO 220 I=1,ND
M=MC(I)
220 SK(M,1)=SK(M,1)*1E8
RETURN
END
SUBROUTINE SOLVE(SK,Q,MS,MB) DIMENSION SK(MS,MB),Q(MS)
K1=MS-1
DO 125 K=1,K1
IF(K+MB-1-MS)105,106,106
105 N=K+MB-1
GOTO 110
106 N=MS
110 I1=K+1
DO 125 I=I1,N
L=I-K+1
C=SK(K,L)/SK(K,1)
J1=MB-L+1
DO 122 J=1,J1
122 SK(I,J)=SK(I,J)-C*SK(K,M)
125 Q(I)=Q(I)-C*Q(K)
Q(MS)=Q(MS)/SK(MS,1)
M=MS-1
DO 145 I1=1,M
I=MS-I1
IF(MS-I+1-MB)135,136,136
135 N=MS-I+1
GOTO 140
136 N=MB
140 DO 142 J=2,N
L=J+I-1
142 Q(I)=Q(I)-SK(I,J)*Q(L)
145 Q(I)=Q(I)/SK(I,1)
WRITE(*,147)
147 FORMAT(5X,'DISPLACEMENT')
WRITE(*,150)(Q(I),I=1,MS)
150 FORMAT(2X,6E14.5)
RETURN
END
SUBROUTINE STRES(Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T) DIMENSION Q(MS),QE(ME),NX(NE,LE),XY(LS,NN),XYE(LS,NE)
DO 400 L=1,LE
DO 160 I=1,NE
DO 160 J=1,NF
N=NF*(I-1)+J
N1=NF*(NX(I,L)-1)+J
160 QE(N)=Q(N1)
WRITE(*,165)L
WRITE(*,170)(QE(I),I=1,ME)
165 FORMAT(4X,'L=',I4)
170 FORMAT(6E14.5)
DO 200 J=1,NE
LJ=NX(J,L)
DO 200 I=1,LS
200 XYE(I,J)=XY(I,LJ)
CALL STE(XYE,QE,NE,LS,ME,E,UM,T)
400 CONTINUE
RETURN
END
SUBROUTINE STIFE(EK,XYE,ME,NE,NF,LS,E,UM,T)
DIMENSION EK(ME,ME),XYE(LS,NE),B(3),C(3)
B(1)=XYE(2,2)-XYE(2,3)
B(2)=XYE(2,3)-XYE(2,1)
B(3)=XYE(2,1)-XYE(2,2)
C(1)=XYE(1,3)-XYE(1,2)
C(2)=XYE(1,1)-XYE(1,3)
C(3)=XYE(1,2)-XYE(1,1)
AE=(B(2)*C(3)-B(3)*C(2))/2
D=E*T/4/(1-UM**2)/AE
DO 30 I=1,3
DO 30 J=1,3
M=2*(I-1)+1
N=2*(J-1)+1
EK(M,N)=D*(B(I)*B(J)+C(I)*C(J)*(1-UM)/2)
EK(M,N+1)=D*(UM*B(I)*C(J)+C(I)*B(J)*(1-UM)/2)
EK(M+1,N)=D*(UM*C(I)*B(J)+B(I)*C(J)*(1-UM)/2)
30 EK(M+1,N+1)=D*(C(I)*C(J)+B(I)*B(J)*(1-UM)/2) RETURN
END
SUBROUTINE STE(XYE,QE,NE,LS,ME,E,UM,T) DIMENSION XYE(LS,NE),QE(ME),SG(3),B(3),C(3),S(3,6) B(1)=XYE(2,2)-XYE(2,3)
B(2)=XYE(2,3)-XYE(2,1)
B(3)=XYE(2,1)-XYE(2,2)
C(1)=XYE(1,3)-XYE(1,2)
C(2)=XYE(1,1)-XYE(1,3)
C(3)=XYE(1,2)-XYE(1,1)
AE=(B(2)*C(3)-B(3)*C(2))/2
D=E*T/2/(1-UM**2)/AE
DO 220 I=1,3
S(1,2*I-1)=D*B(I)
S(2,2*I-1)=D*UM*B(I)
S(3,2*I-1)=D*(1-UM)*C(I)/2
S(1,2*I)=D*UM*C(I)
S(2,2*I)=D*C(I)
220 S(3,2*I)=D*(1-UM)*B(I)/2
DO 250 I=1,3
SG(I)=0
DO 250 K=1,ME
250 SG(I)=SG(I)+S(I,K)*QE(K)
WRITE(*,300)SG
300 FORMAT(5X,'SEGMA',3E20.4)
RETURN
END
悬臂梁计算的输入数据文件如下：
2,3,55,14,10,80,2,
2.1e5,0.3,10,
400,80,400,60,400,40,400,20,400,0,
360,80,360,60,360,40,360,20,360,0,
320,80,320,60,320,40,320,20,320,0,
280,80,280,60,280,40,280,20,280,0,
240,80,240,60,240,40,240,20,240,0,
200,80,200,60,200,40,200,20,200,0,
160,80,160,60,160,40,160,20,160,0,
120,80,120,60,120,40,120,20,120,0,
80,80,80,60,80,40,80,20,80,0,
40,80,40,60,40,40,40,20,40,0,
0,80,0,60,0,40,0,20,0,0,
0,-200,0,0,0,0,0,0,0,0,0,-400,
0,0,0,0,0,0,0,0,0,-400,
0,0,0,0,0,0,0,0,0,-400,
0,0,0,0,0,0,0,0,0,-400,
0,0,0,0,0,0,0,0,0,-400,
0,0,0,0,0,0,0,0,0,-400,
0,0,0,0,0,0,0,0,0,-400,
0,0,0,0,0,0,0,0,0,-400,
0,0,0,0,0,0,0,0,0,-400,
0,0,0,0,0,0,0,0,0,-200,0,0,0,0,0,0,0,0,
1,7,2,2,8,3,3,9,4,4,10,5,1,6,7,2,7,8,3,8,9,4,9,10,
6,12,7,7,13,8,8,14,9,9,15,10,6,11,12,7,12,13,8,13,14, 9,14,15,11,17,12,12,18,13,13,19,14,14,20,15,11,16,17, 12,17,18,13,18,19,14,19,20,16,22,17,17,23,18,18,24,19, 19,25,20,16,21,22,17,22,23,18,23,24,19,24,25,21,27,22, 22,28,23,23,29,24,24,30,25,21,26,27,22,27,28,23,28,29, 24,29,30,26,32,27,27,33,28,28,34,29,29,35,30,26,31,32, 27,32,33,28,33,34,29,34,35,31,37,32,32,38,33,33,39,34, 34,40,35,31,36,37,32,37,38,33,38,39,34,39,40,36,42,37, 37,43,38,38,44,39,39,45,40,36,41,42,37,42,43,38,43,44, 39,44,45,41,47,42,42,48,43,43,49,44,44,50,45,41,46,47, 42,47,48,43,48,49,44,49,50,46,52,47,47,53,48,48,54,49, 49,55,50,46,51,52,47,52,53,48,53,54,49,54,55,
101,102,103,104,105,106,107,108,109,110,
程序输出结果文件如下：
2.58E+02 -1.37E+03 1.30E+02 -5.13E+02 1.05E+02 -2.42E+02 6.41E+01 -1.20E+02 4.37E+01 -7.32E+01 -4.13E+01 -1.35E+03 6.76E+01 -5.32E+02 6.46E+01 -2.54E+02 4.65E+01 -1.28E+02 4.21E+01 -8.12E+01 -6.18E+00 -1.36E+03 8.22E+01 -5.27E+02 6.54E+01 -2.53E+02 4.48E+01 -1.29E+02 4.06E+01 -8.21E+01 2.48E+00 -1.36E+03 8.64E+01 -5.26E+02 6.74E+01 -2.52E+02 4.51E+01 -1.28E+02 4.04E+01 -8.17E+01 4.53E+00 -1.36E+03 8.75E+01 -5.26E+02 6.81E+01 -2.51E+02 4.55E+01 -1.28E+02 4.06E+01 -8.14E+01 4.73E+00 -1.36E+03 8.76E+01 -5.25E+02 6.82E+01 -2.51E+02 4.56E+01 -1.28E+02 4.07E+01 -8.13E+01
3.65E+00 -1.36E+03 8.70E+01 -5.25E+02 6.78E+01 -2.51E+02
4.53E+01 -1.28E+02 4.06E+01 -8.12E+01 -6.27E-01 -1.35E+03 8.39E+01 -
5.24E+02
6.61E+01 -2.51E+02 4.43E+01 -1.27E+02 3.98E+01 -8.10E+01 -1.53E+01 -1.35E+03
7.21E+01 -5.20E+02 5.94E+01 -2.48E+02 4.01E+01 -1.26E+02 3.65E+01 -7.98E+01 -4.41E+01 -1.27E+03 2.95E+01 -4.92E+02 3.53E+01 -2.34E+02 2.63E+01 -1.18E+02 2.21E+01 -7.52E+01 3.01E-06 -
8.84E-06 -6.76E-07 -2.63E-07 -2.43E-07 -7.39E-08 -8.58E-08 -1.60E-08 -7.74E-08 -4.72E-09
程序三：四节点矩形薄板单元程序
如下图为两对边简支、另两对边固支的矩形薄板。

其上受有均布压力q的作用，根据结构与荷载的对称性，可取其1/4为计算对象。

已知E=2.1×105 N/mm2 ，
，厚度h=20mm。

将它划分为4×5的矩形网格，其节点及单元编号如图=0.3
2。

建立如图2所示的离散化计算模型。

计算源程序如下：
PROGRAM MAIN
DIMENSION
SK(90,21),EK(12,12),Q(90),MC(35),XY(2,100),XYE(2,4),QE(12),NX(4,30)
OPEN(7,FILE='INP.DAT')
REWIND 7
READ(7,*)NF,NE,NN,MB,ND,LE,LS
READ(7,*)E,UM,T
10 FORMAT(7I5)
12 FORMAT(6F15.2)
WRITE(*,600)NF,NE,NN,MB,ND,LE,LS,E,UM,T
ME=NE*NF
MS=NN*NF
CALL INPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)
WRITE(*,102)((XY(I,J),I=1,LS),J=1,NN)
102 FORMAT(10X,'XY'/,(2X,6F12.1))
WRITE(*,101)(Q(I),I=1,MS)
101 FORMAT(10X,'Q'/,(2X,3F20.3))
WRITE(*,500)((NX(I,J),I=1,NE),J=1,LE)
WRITE(*,400)(MC(I),I=1,ND)
500 FORMAT(10X,'NX'/,(2X,12I6))
600 FORMAT(10X,'NF NE NN MB ND LE LS E UM T'/7(2X,I4),3(2X,F8.4)) 400 FORMAT(10X,'MC'/,(2X,10I6))
!CALL QS(Q,TL,TH,P,LE,NE,MS,NX)
!WRITE(*,120)(Q(I),I=1,MS)
!120 FORMAT(10X,'Q'/(10X,3F20.3))
CALL STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,& UM,T)
CALL SOLVE(SK,Q,MS,MB)
OPEN(9,FILE='OUT.DAT')
REWIND 9
WRITE(9,200)
WRITE(9,250)(Q(I),I=1,MS)
200 FORMAT(5X,'DISPLACEMENT')
250 FORMAT(2X,6E14.5)
CALL STRES(Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T)
STOP 1000
END
SUBROUTINE INPUT(XY,Q,NX,MC,LS,NN,MS,NE,LE,ND)
DIMENSION XY(LS,NN),Q(MS),NX(NE,LE),MC(ND)
READ(7,*)XY
READ(7,*)Q
READ(7,*)NX
READ(7,*)MC
CLOSE(7)
10 FORMAT(6F11.2)
20 FORMAT(12I5)
RETURN
END
!SUBROUTINE QS(Q,TL,TH,P,LE,NE,MS,NX)
!DIMENSION NX(NE,LE),Q(MS)
!DO I=1,MS
!Q(I)=0.0
!ENDDO
!DO J=1,LE
!I1=3*NX(1,J)-2
!I2=3*NX(2,J)-2
!I3=3*NX(3,J)-2
!I4=3*NX(4,J)-2
!Q(I1)=Q(I1)+P*TL*TH
!Q(I1+1)=Q(I1+1)+P*TL*TH*TH/3
!Q(I1+2)=Q(I1+2)-P*TL*TL*TH/3
!Q(I2)=Q(I2)+P*TL*TH
!Q(I2+1)=Q(I2+1)+P*TL*TH*TH/3
!Q(I2+2)=Q(I2+2)+P*TL*TL*TH/3
!Q(I3)=Q(I3)+P*TL*TH
!Q(I3+1)=Q(I3+1)-P*TL*TH*TH/3
!Q(I3+2)=Q(I3+2)+P*TL*TL*TH/3
!Q(I4)=Q(I4)+P*TL*TH
!Q(I4+1)=Q(I4+1)-P*TL*TH*TH/3
!Q(I4+2)=Q(I4+2)-P*TL*TL*TH/3
!ENDDO
!RETURN
!END
SUBROUTINE
STIFS(SK,EK,Q,NX,XY,XYE,MC,MS,MB,ME,ND,LE,NE,NF,NN,LS,E,& UM,T)
DIMENSION SK(MS,MB),EK(ME,ME),NX(NE,LE),MC(ND),XY(LS,NN),& XYE(LS,NE)
DO 35 I=1,MS
DO 35 J=1,MB
35 SK(I,J)=0
DO 200 L=1,LE
DO 40 J=1,NE
LJ=NX(J,L)
DO 40 I=1,LS
40 XYE(I,J)=XY(I,LJ)
DO 50 I=1,ME
DO 50 J=1,ME
50 EK(I,J)=0.0
CALL STIFE(EK,XYE,ME,NE,LS,E,UM,T)
IF(L.EQ.1)WRITE(*,70)EK
70 FORMAT(10X,'EK'/,(6E14.5))
DO 200 I=1,NE
DO 200 II=1,NF
M=NF*(I-1)+II
M1=NF*(NX(I,L)-1)+II
DO 200 J=1,NE
DO 200 JJ=1,NF
N=NF*(J-1)+JJ
N1=NF*(NX(J,L)-1)+JJ
MN=N1-M1+1
IF(MN)200,200,150
150 SK(M1,MN)=SK(M1,MN)+EK(M,N) 200 CONTINUE
DO 220 I=1,ND
M=MC(I)
220 SK(M,1)=SK(M,1)*1E8
RETURN
END
SUBROUTINE SOLVE(SK,Q,MS,MB) DIMENSION SK(MS,MB),Q(MS)
K1=MS-1
DO 125 K=1,K1
IF(K+MB-1-MS)105,106,106
105 N=K+MB-1
GOTO 110
106 N=MS
110 I1=K+1
DO 125 I=I1,N
L=I-K+1
C=SK(K,L)/SK(K,1)
J1=MB-L+1
DO 122 J=1,J1
122 SK(I,J)=SK(I,J)-C*SK(K,M)
125 Q(I)=Q(I)-C*Q(K)
Q(MS)=Q(MS)/SK(MS,1)
M=MS-1
DO 145 I1=1,M
I=MS-I1
IF(MS-I+1-MB)135,136,136
135 N=MS-I+1
GOTO 140
136 N=MB
140 DO 142 J=2,N
L=J+I-1
142 Q(I)=Q(I)-SK(I,J)*Q(L)
145 Q(I)=Q(I)/SK(I,1)
WRITE(*,147)
147 FORMAT(5X,'DISPLACEMENT') WRITE(*,150)(Q(I),I=1,MS)
150 FORMAT(2X,6E14.5)
RETURN
END
SUBROUTINE STRES(Q,QE,NX,XY,XYE,MS,ME,NE,LE,NF,NN,LS,E,UM,T) DIMENSION Q(MS),QE(ME),NX(NE,LE),XY(LS,NN),XYE(LS,NE)
DO 400 L=1,LE
DO 160 I=1,NE
DO 160 J=1,NF
N=NF*(I-1)+J
N1=NF*(NX(I,L)-1)+J
160 QE(N)=Q(N1)
WRITE(*,165)L
WRITE(*,170)(QE(I),I=1,ME)
165 FORMAT(4X,'L=',I4)
170 FORMAT(6E14.5)
DO 200 J=1,NE
LJ=NX(J,L)
DO 200 I=1,LS
200 XYE(I,J)=XY(I,LJ)
CALL STE(XYE,QE,NE,LS,ME,E,UM,T)
400 CONTINUE
RETURN
END
SUBROUTINE STIFE(EK,XYE,ME,NE,LS,E,UM,T)
DIMENSION EK(ME,ME),XYE(LS,NE)
A=(XYE(1,4)-XYE(1,1))/2
B=(XYE(2,1)-XYE(2,2))/2
A2=A*A
B2=B*B
AB=A2/B2
ET=E*T*T*T/360/(1-UM*UM)/A/B
EK(1,1)=(21.0-6.0*UM+30.0/AB+30.0*AB)*ET
EK(2,1)=(3.0*B+12.0*UM*B+30.0*AB*B)*ET
EK(2,2)=(8.0*B2-8.0*UM*B2+40.0*A2)*ET
EK(3,1)=(-3.0*A-12.0*UM*A-30.0*A/AB)*ET
EK(3,2)=(-30.0*UM*A*B)*ET
EK(3,3)=(8.0*A2-8.0*UM*A2+40.0*B2)*ET
EK(4,1)=(-21.0+6.0*UM-30.0/AB+15.0*AB)*ET
EK(4,2)=(-3.0*B-12.0*UM*B+15.0*AB*B)*ET
EK(4,3)=(3.0*A-3.0*UM*A+30.0/AB*A)*ET
EK(4,4)=EK(1,1)
EK(5,1)=EK(4,2)
EK(5,2)=(-8.0*B2+8.0*UM*B2+20.0*A2)*ET
EK(5,3)=0.0
EK(5,4)=EK(2,1)
EK(5,5)=EK(2,2)
EK(6,1)=-EK(4,3)
EK(6,2)=0.0
EK(6,3)=(-2.0*A2+2.0*UM*A2+20.0*B2)*ET EK(6,4)=-EK(3,1)
EK(6,5)=-EK(3,2)
EK(6,6)=EK(3,3)
EK(7,1)=(21.0-6.0*UM-15.0/AB-15.0*AB)*ET EK(7,2)=(3.0*B-3.0*UM*B-15.0*AB*B)*ET EK(7,3)=(-3.0*A+3.0*UM*A+15.0*A/AB)*ET EK(7,4)=(-21.0+6*UM+15.0/AB-30.0*AB)*ET EK(7,5)=(-3.0*B+3.0*UM*B-30.0*AB*B)*ET EK(7,7)=EK(1,1)
EK(8,1)=-EK(7,2)
EK(8,2)=(2.0*B2-2.0*UM*B2+10.0*A2)*ET EK(8,3)=0.0
EK(8,4)=-EK(7,5)
EK(8,5)=(-2.0*B2+2.0*UM*B2+20.0*A2)*ET EK(8,6)=0.0
EK(8,7)=-EK(2,1)
EK(8,8)=EK(2,2)
EK(9,1)=-EK(7,3)
EK(9,2)=0.0
EK(9,3)=(2.0*A2-2.0*UM*A2+10.0*B2)*ET EK(9,4)=EK(7,6)
EK(9,5)=0.0
EK(9,6)=(-8.0*A2+8.0*UM*A2+20.0*B2)*ET EK(9,7)=EK(6,4)
EK(9,8)=-EK(6,5)
EK(9,9)=EK(6,6)
EK(10,1)=EK(7,4)
EK(10,2)=EK(7,5)
EK(10,3)=-EK(7,6)
EK(10,4)=EK(7,1)
EK(10,5)=EK(7,2)
EK(10,6)=-EK(7,3)
EK(10,7)=EK(4,1)
EK(10,8)=-EK(4,2)
EK(10,9)=-EK(4,3)
EK(10,10)=EK(1,1)
EK(11,1)=EK(18,4)
EK(11,2)=EK(8,5)
EK(11,3)=0.0
EK(11,4)=EK(8,1)。