有限元实验报告.doc

第一个区域剖分实验及处理结果如下：（1）节点总数ND=15（2）单元总数NE=16（3）单元节点编号数组I3（3，NE）;单元序号12345678910111213141516ijm142524532563547578586896117108711811128129101314111014111415111512（4）节点坐标数组XY（2，ND）节点序号123456789101112131415Xy00010210.511.2122121.52230.731.33240.541.242（5）第一类边界节点数ND1=12（6）第一类边界节点号数组NB1(ND1)和场值数组U1(ND1)第一类 边界节点序号 123456789101112第一类边界节点号场值112.53.334.466.37.259.2310.1412.1713.0814.1115.12给出有限元法程序框图及区域剖分图（如下）：xy1234012345698 710111212151413（1）（2）（3）（4）（5）（6）（7）（8）（9）（10）（11）（13）（14）（15）（16）（12）源程序如下： PROGRAM youxianyuan01 character *20 filename1 PARAMETER(ND=15,NE=16,ND1=12) DIMENSION I3(3,16),XY(2,15),U1(12),NB1(12),SK(15,15),U(15) OPEN(11,FILE=NE.txt,STATUS=old) DO 100 J=1,NE100 READ(11,*)(I3(I,J),I=1,3) READ(11,*)(U1(i),i=1,ND1) DO 2 J=1,ND200 READ(11,*)(XY(I,J),I=1,2) READ(11,*)(NB1(I),I=1,ND1) CLOSE(11) CALL MBW(NE,I3,IW) CALL UK1(ND,NE,IW,I3,XY,SK) CALL UB1(ND1,NB1,U1,ND,IW,SK,U) CALL LDLT(SK,ND,IW,U,IE) write(*,*)请输入数据文件名: read(*,*)filename1 open(6,file=hebin.txt,status=unknown) write(6,510)U close(6)510 format(/3e15.6) END SUBROUTINE MBW(NE,I3,IW) DIMENSION I3(3,NE) IW=0 DO 10 I=1,NE M=MAX(IABS(I3(1,I)-I3(2,I),IABS(I3(2,I)-I3(3,I),& & IABS(I3(3,I)-I3(1,I) IF(M+1.GT.IW)IW=M+110 CONTINUE RETURN END SUBROUTINE UK1(ND,NE,IW,I3,XY,SK) DIMENSION I3(3,NE),XY(2,ND),SK(ND,IW) DIMENSION X(3),Y(3) REAL KE(3,3) DO 10 I=1,ND DO 10 J=1,IW10 SK(I,J)=0. DO 20 L=1,NE DO 30 J=1,3 I=I3(J,L) X(J)=XY(1,I)30 Y(J)=XY(2,I) CALL UKE1(X,Y,KE) DO 40 J=1,3 NJ=I3(J,L) DO 40 K=1,J NK=I3(K,L) IF(NJ.LT.NK)GOTO 50 NK=NK-NJ+IW SK(NJ,NK)=SK(NJ,NK)+KE(J,K) GOTO 40 50 NJ=NJ-NK+IW SK(NK,NJ)=SK(NK,NJ)+KE(J,K) NJ=NJ+NK-IW40 CONTINUE20 CONTINUE RETURN END SUBROUTINE UKE1(X,Y,KE) DIMENSION X(3),Y(3),A(3),B(3) REAL KE(3,3) A(1)=Y(2)-Y(3) A(2)=Y(3)-Y(1) A(3)=Y(1)-Y(2) B(1)=X(3)-X(2) B(2)=X(1)-X(3) B(3)=X(2)-X(1) S=2.*(A(1)*B(2)-A(2)*B(1) DO 10 I=1,3 DO 10 J=1,I10 KE(I,J)=(A(I)*A(J)+B(I)*B(J)/S RETURN END SUBROUTINE UB1(ND1,NB1,U1,ND,IW,SK,U) DIMENSION NB1(ND1),U1(ND1),SK(ND,IW),U(ND) DO 10 I=1,ND10 U(I)=0. DO 20 I=1,ND1 J=NB1(I) SK(J,IW)=SK(J,IW)*1.E1020 U(J)=SK(J,IW)*U1(I) RETURN END SUBROUTINE LDLT(A,N,IW,P,IE) DIMENSION A(N,IW),P(N) DO 15 I=1,N IF(I.LE.IW) GOTO 20 IT=I-IW+1 GOTO 3020 IT=130 K=I-1 IF(I.EQ.1)GOTO 40 DO 25 L=IT,K IL=L+IW-I B=A(I,IL) A(I,IL)=B/A(L,IW) P(I)=P(I)-A(I,IL)*P(L) MI=L+1 DO 25 J=MI,I IJ=J+IW-I JL=L+IW-J25 A(I,IJ)=A(I,IJ)-A(J,JL)*B40 IF(A(I,IW).EQ.0.)GOTO 10015 CONTINUE DO 45 J=1,N IF(J.LE.IW)GOTO 60 IT=N-J+IW GOTO 7060 IT=N70 I=N-J+1 P(I)=P(I)/A(I,IW) IF(J.EQ.1)GOTO 45 K=I+1 DO 65 MJ=K,IT IJ=I-MJ+IW65 P(I)=P(I)-P(MJ)*A(MJ,IJ)45 CONTINUE IE=0 GOTO 110100 IE=1110 RETURN END运行时输入的数据：1,4,25,2,45,3,25,6,35,4,75,7,85,8,68,9,6,11,7,108,7,118,11,128,12,910,13,1411,10,1411,14,1511,15,121,0.5,0.33,0.46,0.3,0.25,0.23,0.14,0.17,0.08,0.11,0.120,00,10,21,0.51,1.21,22,12,1.52,23,0.73,1.33,24,0.54,1.24,21,2,3,4,6,7,9,10,12,13,14,15输出结果： 0.100000E+01 0.500000E+00 0.330000E+00 0.460000E+00 0.375833E+00 0.300000E+00 0.250000E+00 0.243475E+00 0.230000E+00 0.140000E+00 0.161536E+00 0.170000E+00 0.800000E-01 0.110000E+00 0.120000E+00第二个区域剖分实验及处理结果如下：（1） 节点总数ND=21（2） 单元总数NE=28（3） 单元节点编号数组I3(3,NE)单元序号12345678910111213141516171819202122232425262728IJm1451527542576273268673687412117129711879911108910121314111214151114151416101115101516132120192021141320142019141918171819161418161817（4） 节点坐标数组XY(2,ND)节点序号123456789101112131415161718192021Xy00010211.50.50.750.51.511.2121.51.752221.52130.731.32.51.7532423.51.6541.23.50.9540.5（5） 第一类边界节点总数ND1=12（6） 第一类边界节点号数组NB1(ND1)和场值数组U1(ND1)第一类边界点序号123456789101112第一类边界节点号场值1121.530.3340.4680.3100.25120.23130.14160，17170.08190.11210.12给出有限元法程序框图及区域剖分图（如下）：xy12340123456789101112131415161718192021（1）（2）（3）（4）（5）（6）（7）（8）（9）（10）（11）（12）（13）（14）（15）（16）（17）（18）（20）（21）（22）（23）（24）（25）（26）（27）（28）（19）源程序如下： PROGRAM youxianyuan02 character *20 filename1 PARAMETER(ND=21,NE=28,ND1=12) DIMENSION I3(3,28),XY(2,21),U1(12),NB1(12),SK(21,21),U(21) OPEN(11,FILE=NE.txt,STATUS=old) DO 100 J=1,NE100 READ(11,*)(I3(I,J),I=1,3) READ(11,*)(U1(i),i=1,ND1) DO 200 J=1,ND200 READ(11,*)(XY(I,J),I=1,2) READ(11,*)(NB1(I),I=1,ND1) CLOSE(11) CALL MBW(NE,I3,IW) CALL UK1(ND,NE,IW,I3,XY,SK) CALL UB1(ND1,NB1,U1,ND,IW,SK,U) CALL LDLT(SK,ND,IW,U,IE) write(*,*)请输入数据文件名: read(*,*)filename1 open(6,file=hebin.txt,status=unknown) write(6,510)U close(6)510 format(/3e15.6) END SUBROUTINE MBW(NE,I3,IW) DIMENSION I3(3,NE) IW=0 DO 10 I=1,NE M=MAX(IABS(I3(1,I)-I3(2,I),IABS(I3(2,I)-I3(3,I),& & IABS(I3(3,I)-I3(1,I) IF(M+1.GT.IW)IW=M+110 CONTINUE RETURN END SUBROUTINE UK1(ND,NE,IW,I3,XY,SK) DIMENSION I3(3,NE),XY(2,ND),SK(ND,IW) DIMENSION X(3),Y(3) REAL KE(3,3) DO 10 I=1,ND DO 10 J=1,IW10 SK(I,J)=0. DO 20 L=1,NE DO 30 J=1,3 I=I3(J,L) X(J)=XY(1,I)30 Y(J)=XY(2,I) CALL UKE1(X,Y,KE) DO 40 J=1,3 NJ=I3(J,L) DO 40 K=1,J NK=I3(K,L) IF(NJ.LT.NK)GOTO 50 NK=NK-NJ+IW SK(NJ,NK)=SK(NJ,NK)+KE(J,K) GOTO 40 50 NJ=NJ-NK+IW SK(NK,NJ)=SK(NK,NJ)+KE(J,K) NJ=NJ+NK-IW40 CONTINUE20 CONTINUE RETURN END SUBROUTINE UKE1(X,Y,KE) DIMENSION X(3),Y(3),A(3),B(3) REAL KE(3,3) A(1)=Y(2)-Y(3) A(2)=Y(3)-Y(1) A(3)=Y(1)-Y(2) B(1)=X(3)-X(2) B(2)=X(1)-X(3) B(3)=X(2)-X(1) S=2.*(A(1)*B(2)-A(2)*B(1) DO 10 I=1,3 DO 10 J=1,I10 KE(I,J)=(A(I)*A(J)+B(I)*B(J)/S RETURN END SUBROUTINE UB1(ND1,NB1,U1,ND,IW,SK,U) DIMENSION NB1(ND1),U1(ND1),SK(ND,IW),U(ND) DO 10 I=1,ND10 U(I)=0. DO 20 I=1,ND1 J=NB1(I) SK(J,IW)=SK(J,IW)*1.E1020 U(J)=SK(J,IW)*U1(I) RETURN END SUBROUTINE LDLT(A,N,IW,P,IE) DIMENSION A(N,IW),P(N) DO 15 I=1,N IF(I.LE.IW) GOTO 20 IT=I-IW+1 GOTO 3020 IT=130 K=I-1 IF(I.EQ.1)GOTO 40 DO 25 L=IT,K IL=L+IW-I B=A(I,IL) A(I,IL)=B/A(L,IW) P(I)=P(I)-A(I,IL)*P(L) MI=L+1 DO 25 J=MI,I IJ=J+IW-I JL=L+IW-J25 A(I,IJ)=A(I,IJ)-A(J,JL)*B40 IF(A(I,IW).EQ.0.)GOTO 10015 CONTINUE DO 45 J=1,N IF(J.LE.IW)GOTO 60 IT=N-J+IW GOTO 7060 IT=N70 I=N-J+1 P(I)=P(I)/A(I,IW) IF(J.EQ.1)GOTO 45 K=I+1 DO 65 MJ=K,IT IJ=I-MJ+IW65 P(I)=P(I)-P(MJ)*A(MJ,IJ)45 CONTINUE IE=0 GOTO 110100 IE=1110 RETURN END输入数据如下：1 4 51 5 2 4 7 55 7 27 6 26 3 26 7 86 8 34 12 77 12 117 11 97 9 89 11 109 10 812 13 1412 14 1111 14 1515 14 1615 10 1115 16 1013 21 20 20 21 1913