random.for

随机有限元分析的经典程序

          CCTT=C00-C01+C02-C03-C10+C11-C12+C13+C20
               -C21+C22-C23-C30+C31-C32+C33
          COV(M,N)=CCTT*SQRT(VAR(LV1)*VAR
      ﹡  (LV2))/TX1/TY1/TX2/TY2/4.0
160       CONTINUE
          CONTINUE
          RETURN
          END

          SUBROUTINE UDIAVAR(I,J,LM,LN,LV1,LV2,SCB)
          COMMON/FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
          COMMON/CONTR1/NEAD(2), NEADA(2), NEA, NEB, MY(2,2),NEL
          COMMON /TRCOV/COV(150,150), DIB3(150),
      ﹡  PSI(150,150), DIA3(150)
          COMMON/TRCOV1/VAR(4), SCA(10), 
      ﹡  INEA(150,5), AERA(150)
          COMMON/TRCOV2/P(4), H(4), SP(4), DSP(2,4)
          B1=100000.0
          B2=B1
          DO 30 K1=1,4
          L1=INEA(LM,K1)
          DO 30 K2=1,4
          L2=INEA(LN,K2)
          B3=ABS(X(L1)-X(L2))
          B4=ABS(Y(L1)-Y(L2))
          IF(B3.LT.B1) B1=B3
          IF(B4.LT.B2) B2=B4
30        CONTINUE
          IF(B1.GT.SCB.OR.B2.GT.SCB) GOTO 90
          DO 50 M1=1,4
          DO 50 M2=1,4
          CALL SHAPE(M1,M2)
          CALL JACOB(LM,XJI)
          CALL SF(LM,U1,U2)
          DO 50 M3=1,4
          DO 50 M4=1,4
          CALL SHAPE(M3,M4)
          CALL TACOB(LN,XJJ)
          CALL SF(LN,U3,U4)
          VA=SCB-ABS(U1-U3)
          VB=SCB-ABS(U2-U4)
          IF(VA.LT.0.0.OR.VB.LT.0.0) GOTO 50
          COV(I,J)=COV(I,J)+VA*VB*H(M1)*H(M2)
      ﹡  *H(M3)*H(M4)*XJI*XJJ*SQRT
      ﹡  (VAR(LV1)*VAR(LV2))/AERA(LN)/AERA(LM)/SCB/SCB
50        CONTINUE
90        RETURN
          END
 
          SUBROUTINE SF(M,UA1,UA2)
          COMMON/FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
          COMMON/TRCOV1/VAR(4), SCA(10), 
      ﹡  INEA(150,5), AERA(150)
          COMMON/TRCOV2/P(4), H(4), SP(4), DSP(2,4)
          UA1=0.0
          UA2=0.0
          DO 10 K=1,4
          L=INEA(M,K)
          UA1=UA1+SP(K)*X(L)
10        UA2=UA2+SP(K)*Y(L)
          RETURN
          END
  
          SUBROUTINE SHAPE(M,N)
          COMMON/FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
          COMMON/TRCOV1/VAR(4), SCA(10), 
      ﹡  INEA(150,5), AERA(150)
          COMMON/TRCOV2/P(4), H(4), SP(4), DSP(2,4)
          SP(1)=(1-P(M)-P(N)+P(M)*P(N))/4
          SP(2)=(1+P(M)-P(N)-P(M)*P(N))/4
          SP(3)=(1+P(M)+P(N)+P(M)*P(N))/4
          SP(4)=(1-P(M)+P(N)-P(M)*P(N))/4
          DSP(1,1)=(-1+P(N))/4
          DSP(1,2)=(1-P(N))/4
          DSP(1,3)=(1+P(N))/4
          DSP(1,4)=(-1-P(N))/4
          DSP(2,1)=(-1+P(M))/4
          DSP(2,2)=(-1-P(M))/4
          DSP(2,3)=(1+P(M))/4
          DSP(2,4)=(1-P(M))/4
          RETURN
          END
   
          SUBROUTINE JACOB(I,XJI)
          DIMENSION UB(4)
          COMMON/FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
          COMMON/TRCOV1/VAR(4), SCA(10), 
      ﹡  INEA(150,5), AERA(150)
          COMMON/TRCOV2/P(4), H(4), SP(4), DSP(2,4)
          DO 10 K=1,4
10        UB(K)=0.0
          DO 20 K=1,4
          L=INEA(I,K)
          UB(1)=UB(1)+DSP(1,K)*X(L)
          UB(2)=UB(2)+DSP(2,K)*Y(L)
          UB(3)=UB(3)+DSP(2,K)*X(L)
20        UB(4)=UB(4)+DSP(1,K)*Y(L)
          XJI=UB(1)*UB(2)-UB(3)*UB(4)
          RETURN
          END
 
          SUBROUTINE DVM(I,LV,LE)
          COMMON/CONTR1/NEAD(2), NEADA(2), NEA, NEB, MY(2,2),NEL
          IF(I,GT.NEADA(1)) GOTO 20
          IF(I.GT.NEAD(1)) GOTO 10
          LV=1
          LE=I
          IF(MY(1,1).EQ.0) LV=2  
      ﹡  GOTO 40
10        LV=2
          LE=I-NEAD(1)
          GOTO 40
20        CONTINUE
          IF(I.GT.(NEADA(1)+NEAD(2))) GOTO 30
          LV=3
          LE=I-(MY(1,1)+MY(1,2)-1)*NEAD(1)
          IF(MY(2,1).EQ.0) LV=4
          GOTO 40
30        LV=4
          LE=I-(MY(1,1)+MY(1,2)-1)*NEAD(1)-NEAD(2)
40        CONTINUE
          RETURN
          END
 
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
                  APPROXIMATE TRANSFORMATION
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE ERLE
          COMMON /TRCOV/COV(150,150), DIB3(150), 
      ﹡  PSI(150,150), DIA3(150)
      ﹡  /CONTR1/NEAD(2), NEADA(2), NEA, NEB, MY(2,2), NEL
          A=0.0
          NNP=NEA+NEL
          DO 10 J=1,NNP
          DO 10 J=1,NNP
10        A=A+COV(I,J)*COV(I,J)
          A=SQRT(A)
          S=FLOAT(NNP)/4.0
          NEB=IFIX(S)
20        CONTINUE
          B=0.0
          DO 30 I=1,NNP
          DO 30 J=1,NNP
          DO 30 K=1,NEB
30        B=B+(PSI(I,K)*DIA3(K)*PSI(J,K)**2
          B=SQRT(B)
          C=ABS((A-B)/A)*100.0
          IF(C.LT.5.) GOTO 40
          NES=NEB+2
          GOTO 20
40        WRITE(16,50) NEB,C
50        FORMAT(/5X,'NEB=',I3,10X,'ERROR=',F10.7,'%'/)
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
     THE SECOND ORDER STATISTICS OF THE NODAL DISPLACEMENT
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE DISPL
          COMMON /CONTR/NX, NE, MT, NR, MP, MTX,
      ﹡  NN, LNY, NERL, NYC
      ﹡  /CONTR1/NEAD(2), NEADA(2), NEA, NEB, MY(2,2), NEL
      ﹡  /SUIJI1/VB(600), VC(600,150), V2(600), F(600), QF(600)
      ﹡  /TRCOV/COV(150,150), DIB3(150), PSI(150,150), DIA3(150)
          WRITE(16,10)
10        FORMAT(/1X,9('﹡').'THE MEANS AND VARIANCES 
      ﹡  OF DISPLACEMENT MENT IN EVERY NODE',9('﹡')/)
          DO 200 I=1,NX
          DO 200 K=1,2
          II=(I-1)*2+K
          ED=VB(II)+V2(II)
          VD=0.0
          DO 100 J=1,NEB
          VD=VD+VC(II,J)**2*DIA3(J)]
100       CONTINUE
          IF(K.EQ.1) WRITE(16,88)I,ED,VD
          IF(K.EQ.2) WRITE(16,99)ED,VD
88        FORMAT(1X,'NO.',I5,8X,'X:',8X,'EV=',
      ﹡  E12.5,8X,'VAR V=',E12.5)
99        FORMAT(17X,'Y:',8X,'EV=',E12.5,8X,
      ﹡  'VAR V=',E12.5)
200       CONTINUE
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          THE SECOND ORDER STATISTICS OF ELEMENT STRESS
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE STOCH
          DIMENSION UV1(2,6),UV2(2,6),UVC(2,6,150),PVD(6,6)
          COMMON /CONTR/NX, NE, MT, NR, MP, MTX,
      ﹡  NN, LNY, NERL, NYC
      ﹡  /CONTR1/NEAD(2), NEADA(2), NEA, NEB, MY(2,2), NEL
      ﹡  /FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
      ﹡  /SUIJI1/VB(600), VC(600,150), V2(600), F(600), QF(600)
      ﹡  /SUIJI2/D(3,3,2), DD(3,3,4), DD2(3,3,4), 
      ﹡  DD2EG(3,3,2), EK(6,6)
      ﹡  /SUIJI3/BB(3,6,400), AREA(400)
      ﹡  /SUIJI4/MDA(600), MXB, MN
          /TRCOV/COV(150,150), DIB3(150), 
      ﹡  PSI(150,150), DIA3(150)
      ﹡  /SMCOV/SMEAN(3,400), COVAR(3,3,400)
          DO 800 IE=1,NE
          CALL DSL(IE,MTY,INT,INS,MT1,MT2)
          II=INE(IE,1)
          JJ=INE(IE,2)
          KK=INE(IE,3)
          UV1(1,1)=VB(2*II-1)
          UV1(1,2)=VB(2*II)
          UV1(1,3)=VB(2*JJ-1)
          UV1(1,4)=VB(2*JJ)
          UV1(1,5)=VB(2*KK-1)
          UV1(1,6)=VB(2*KK)
          UV2(1,1)=VB(2*II-1)
          UV2(1,2)=VB(2*II)
          UV2(1,3)=VB(2*JJ-1)
          UV2(1,4)=VB(2*JJ)
          UV2(1,5)=VB(2*KK-1)
          UV2(1,6)=VB(2*KK)
          DO 100 M=1,3
          SMEAN(M,IE)=0.0
          DO 111 K=1,6
111       PVD(1,K)=0.0
          DO 100 N=1,3
          DO 110 L=1,6
110       PVD(1,N)=PVD(1,N)+BB(N,L,IE)*(UV1(1,L)+UV2(1,L))
100       SMEAN(M,IE)=SMEAN(M,IE)+D(M,N,MTY)*PVD(1,N)/AREA(IE)
          DO 112 K=1,6
          DO 112 L=1,6
112       PVD(K,L)=0.0
          DO 120 J=1,NEB
          UV1(1,1)=VB(2*II-1)
          UV1(1,2)=VB(2*II)
          UV1(1,3)=VB(2*JJ-1)
          UV1(1,4)=VB(2*JJ)
          UV1(1,5)=VB(2*KK-1)
          UV1(1,6)=VB(2*KK)
          UVC(1,1,J)=VB(2*II-1,J)
          UVC(1,2,J)=VB(2*II,J)
          UVC(1,3,J)=VB(2*JJ-1,J)
          UVC(1,4,J)=VB(2*JJ,J)
          UVC(1,5,J)=VB(2*KK-1,J)
          UVC(1,6,J)=VB(2*KK,J)
120       CONTINUE
          DO 130 M=1,6
          DO 130 J=1,NEB
          PVD(1,M)=PVD(1,M)+UVC(1,M,J)*PSI(INT,J)*DIA3(J)
          PVD(2,M)=PVD(2,M)+0.5*UV1(1,M)*PSI(INT,J)*PSI(INT,J)*DIA3(J)
          IF(INS.EQ.0) GOTO 130
          PVD(3,M)=PVD(3,M)+UVC(1,M,J)*PSI(INS,J)*DIA3(J)
          PVD(4,M)=PVD(4,M)+0.5*UV1(1,M)*PSI(INS,J)*PSI(INS,J)*DIA3(J)
          PVD(5,M)=PVD(5,M)+UV1(1,M)*PSI(INT,J)*PSI(INS,J)*DIA3(J)
130       CONTINUE
          DO 140 M=1,3
          DO 140 N=1,3
          DO 140 L=1,6
          RR1=DD(M,N,MT1)*BB(N,L,IE)*PVD(1,L)
          RR2=DD(M,N,MT2)*BB(N,L,IE)*PVD(3,L)
          RR3=DD2(M,N,MT1)*BB(N,L,IE)*PVD(2,L)
          RR4=DD2(M,N,MT2)*BB(N,L,IE)*PVD(4,L)
          RR5=DD2EG(M,N,MTY)*BB(N,L,IE)*PVD(5,L)
          SMEAN(M,IE)=SMEAN(M,IE)+(RR1+RR2+RR3+RR4+RR5)/AREA(IE)
140       CONTINUE
          WRITE(16,400) IE,(SMEAN(I,IE),I=1,3)
400       FORMAT(30X,'************',I3,'************'/5X,
      ﹡  'MEANX=',E12.5,4X,'MEANY=',E12.5,4X,'MEANXY=',E12.5/)
          JE=IE
          CALL DSL(JE,NTY,JNT,JNS,NT1,NT2)
          II=INE(JE,1)
          JJ=INE(JE,2)
          KK=INE(JE,3)
          UV1(2,1)=VB(2*II-1)
          UV1(2,2)=VB(2*II)
          UV1(2,3)=VB(2*JJ-1)
          UV1(2,4)=VB(2*JJ)
          UV1(2,5)=VB(2*KK-1)
          UV1(2,6)=VB(2*KK)
          DO 150 J=1,NEB
          UVC(2,1,J)=VB(2*II-1,J)
          UVC(2,2,J)=VB(2*II,J)
          UVC(2,3,J)=VB(2*JJ-1,J)
          UVC(2,4,J)=VB(2*JJ,J)
          UVC(2,5,J)=VB(2*KK-1,J)
          UVC(2,6,J)=VB(2*KK,J)
150       CONTINUE
          DO 153 KM=1,3
          DO 153 KN=1,3
153       COVAR(KM,KN,JE)=0.0
          DO 200 J=1,NEB
          DO 155 K=1,6
          DO 155 L=1,6
155       PVD(K,L)=0.0
          DO 160 M=1,3
          DO 160 N=1,3
          DO 160 L=1,6
          PVD(1,M)=PVD(1,M)+DD(M,N,MT1)*BB(N,L,IE)
      ﹡  *UV1(1,L)*PSI(INT,J)/AREA(IE) 
          PVD(2,M)=PVD(2,M)+DD(M,N,MT2)*BB(N,L,IE)
      ﹡  *UV1(1,L)*PSI(INS,J)/AREA(IE)
          PVD(3,M)=PVD(3,M)+D(M,N,MTY)*BB(N,L,IE)
      ﹡  *UVC(1,L,J)/AREA(IE)
          PVD(4,M)=PVD(4,M)+DD(M,N,NT1)*BB(N,L,JE)
      ﹡  *UV1(2,L)*PSI(JNT,J)/AREA(JE)
          PVD(5,M)=PVD(5,M)+DD(M,N,NT2)*BB(N,L,JE)
      ﹡  *UV1(2,L)*PSI(JNS,J)/AREA(JE)
          PVD(6,M)=PVD(6,M)+D(M,N,MTY)*BB(N,L,JE)
      ﹡  *UVC(2,L,J)/AREA(JE)
160       CONTINUE
          DO 170 M=1,3
          DO 170 N=M,3
          RR1=PVD(1,M)*PVD(4,N)+PVD(1,M)*PVD(5,N)
          RR2=PVD(1,M)*PVD(6,N)+PVD(2,M)*PVD(4,N)
          RR3=PVD(2,M)*PVD(5,N)+PVD(2,M)*PVD(6,N)
          RR4=PVD(3,M)*PVD(4,N)+PVD(3,M)*PVD(5,N)
          RR5=PVD(3,M)*PVD(6,N)
170       COVAR(M,N,JE)=COVAR(M,N,JE)+
      ﹡  (RR1+RR2+RR3+RR4+RR5)*DIA3(J)
200       CONTINUE 
          WRITE(16,500) ((COVAR(I,J,JE),J=1,3),I=1,3)
500       FORMAT(5X,'XX^XX=',E12.5,4X,'XX^YY='.E12.5,4X,
      ﹡  'XX^YY=',E12.5/27X,'YY^YY=',E12.5,4X,'YY^XY=',
      ﹡  E12.5/49X,'XY^XY=',E12.5/)
800       CONTINUE 
          RETURN
          END

          SUBROUTINE DSL(IE,MTY,INT,INS,MT1,MT2)
          COMMON /CONTR/NX, NE, MT, NR, MP, MTX,
      ﹡  NN, LNY, NERL, NYC
      ﹡  /CONTR1/NEAD(2), NEADA(2), NEA, NEB, MY(2,2), NEL
          MTY=INE(IE,4)
          INT=INE(IE,5)+(MTY-1)*NEADA(1)
          INS=(MY(MTY,1)+MY(MTY,2)-1)*(INT+NEAD(MTY))
          MT1=2*(MTY-1)+1
          MT2=MT1+1
          IF(INS.EQ.0.AND.MY(MTY,1).EQ.0)
      ﹡  INS=INT
          RETURN
          END
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
      THE MEANS AND VARIANCES OF STRESS IN EVERY NODE
﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡﹡
          SUBROUTINE SETOUT
          DIMENSION XX(3)
          COMMON /CONTR/NX, NE, MT, NR, MP, MTX,
      ﹡  NN, LNY, NERL, NYC
          COMMON/SMCOV/SMEAN(3,400), COVAR(3,3,400)
      ﹡  /FONDA/AE(2,8), X(300), Y(300), INE(400,5), 
      ﹡  RR(40), NPP(30), PP(60)
      ﹡  /SUIJI4/MDA(600), MXB, MN
          WRITE(16,10)
10        FORMAT(15X,'THE MEANS AND VARIANCES OF 
      ﹡  STRESS IN EVERY NODE'/)
          WRITE(16,20)
20        FORMAT(1X/5X,'NO.',10X,'E(SIGMAX)',12X,
      ﹡  'E(SIGMAY)',12X,'E(SIGMAXY)'/)
          DO 40 J=1,2
          IF(J.EQ.1) GOTO 21
          WRITE(16,22)
22        FORMAT(1X/5X,'NO.',9X,'VAR(SIGMAX)',9X,
      ﹡  'VAR(SIGMAY)',9X,'VAR(SIGMAXY)'/)
21        CONTINUE
          DO 40 I=1,NX
          DO 30 MK=1,3
30        XX(MK)=0.0
          M=0
          DO 90 IK=1,NE
          DO 70 IL=1,3
          IF(ABS(I-INE(IK,IL)).GT.MXB) GOTO 90
          IF(INE(IK,IL).EQ.I) GOTO 80
70        CONTINUE
          GOTO 90
80        M=M+1
          DO 100 IM=1,3
          IF(J.EQ.2) GOTO 110
          XX(IM)=XX(IM)+SMEAN(IM,IK)
          GOTO 100
110       XX(IM)=XX(IM)+COVAR(IM,IM,IK)
100       CONTINUE
90        CONTINUE
          DO 120 IM=1,3
120       XX(IM)=XX(IM)/M
          WRITE(16,140) I,(XX(IM),IM=1,3)
140       FORMAT(5X,I3,8X,F13.6,8X,F13.6,8X,F13.6)
40        CONTINUE
          RETURN
          END