apd_2.for

fortran写的多层弹性体系的弹性力学求解

        FV(1)=V6*P1
        FV(2)=V6*P2
        FV(3)=V7*UU
        FV(4)=V6*P3
        FV(6)=V9*UU
        FV(5)=V4*P4
        FV(7)=V8*WW
        DO 80 K=1,7
80      YV(ISW,K)=YV(ISW,K)+FV(K)*WR
        IF(Z.EQ.0) THEN
        SV(2)=V7
        YW(ISW,2)=YW(ISW,2)+SV(2)*WR
        ENDIF
        ENDIF
        IF(H.EQ.1) THEN
        FH(1)=V4*PPR
        FH(2)=V4*PQ
        FH(3)=V4*PZ
        FH(4)=V4*PT
        FH(5)=0
        IF(RD.NE.0) FH(5)=V5*U2/Q
        FH(6)=V5*H2
        FH(7)=V6*H0
        FH(8)=V8*U0
        FH(9)=V9*PW
        FH(10)=V3*U2
        DO 90 K=1,10
90      YH(ISW,K)=YH(ISW,K)+FH(K)*WR
        IF(Z.EQ.0) THEN
        SH(1)=V4
        SH(3)=V3
        SH(2)=0
        IF(RD.NE.0) SH(2)=V5/Q
        DO 100 K=1,3
100     YI(ISW,K)=YI(ISW,K)+SH(K)*WR
        ENDIF
        ENDIF
30      CONTINUE
        DO 40 ISW=1,2
        IF(Z.EQ.0) THEN
        DO 110 K=3,4
110     YW(ISW,K)=YW(ISW,K)*C
        YW(ISW,1)=YW(ISW,1)*C
        ENDIF
        IF(V.EQ.1) THEN
        DO 120 K=1,7
120     YV(ISW,K)=YV(ISW,K)*C
        IF(Z.EQ.0)YW(ISW,2)=YW(ISW,2)*C
        CR(ISW)=YV(ISW,1)-YV(ISW,3)
        CQ(ISW)=YV(ISW,2)*R1+YV(ISW,3)
        CZ(ISW)=-YV(ISW,4)
        TR(ISW)=0
        TQ(ISW)=0
        TZ(ISW)=YV(ISW,5)
        UV(ISW)=.5*R5*YV(ISW,6)
        VV(ISW)=0
        WV(ISW)=-.5*R5*YV(ISW,7)
        ENDIF
        IF(H.EQ.1) THEN
        DO 130 K=1,10
        YH(ISW,K)=YH(ISW,K)*C
        IF((K.LE.3).AND.(Z.EQ.0)) YI(ISW,K)=YI(ISW,K)*C
130     CONTINUE
        CA(ISW)=YH(ISW,1)-YH(ISW,5)
        CB(ISW)=R1*YH(ISW,2)+YH(ISW,5)
        CD(ISW)=-YH(ISW,3)
        TD(ISW)=YH(ISW,4)-YH(ISW,5)
        TE(ISW)=.5*(YH(ISW,6)+YH(ISW,7))
        TF(ISW)=.5*(YH(ISW,6)-YH(ISW,7))
        UH(ISW)=.25*R5*(YH(ISW,10)-YH(ISW,8))
        VH(ISW)=.25*R5*(YH(ISW,10)+YH(ISW,8))
        WH(ISW)=-.5*R5*YH(ISW,9)
        ENDIF
        IF(Z.EQ.0) THEN
        CALL SURFA(ISW,PO,RR,YW,YI,H,V,SSU,SSV)
        IF(V.EQ.1) THEN
        CR(ISW)=CR(ISW)+SSU(1)
        CQ(ISW)=CQ(ISW)+SSU(2)
        CZ(ISW)=CZ(ISW)+SSU(3)
        UV(ISW)=UV(ISW)+SSU(4)
        WV(ISW)=WV(ISW)+SSU(5)
        ENDIF
        IF(H.EQ.1) THEN
        CA(ISW)=CA(ISW)+SSV(1)
        CB(ISW)=CB(ISW)+SSV(2)
        TD(ISW)=TD(ISW)+SSV(3)
        TE(ISW)=TE(ISW)+SSV(4)
        TF(ISW)=TF(ISW)+SSV(5)
        UH(ISW)=UH(ISW)+SSV(6)
        VH(ISW)=VH(ISW)+SSV(7)
        WH(ISW)=WH(ISW)+SSV(8)
        ENDIF
        ENDIF
        IF(V.EQ.1) THEN
        ER(ISW)=CR(ISW)-(CQ(ISW)+CZ(ISW))*PO
        EQ(ISW)=CQ(ISW)-(CZ(ISW)+CR(ISW))*PO
        EZ(ISW)=CZ(ISW)-(CR(ISW)+CQ(ISW))*PO
        ENDIF
        IF(H.EQ.1) THEN
        EG(ISW)=CA(ISW)-(CB(ISW)+CD(ISW))*PO
        EH(ISW)=CB(ISW)-(CD(ISW)+CA(ISW))*PO
        EI(ISW)=CD(ISW)-(CA(ISW)+CB(ISW))*PO
        ENDIF
40      CONTINUE
        RETURN
        END
*****************************************************
        SUBROUTINE DOUBLE(V,H,F,R1,CE,D)
        COMMON /CCC5/CR(2),CQ(2),CZ(2),TR(2),TQ(2),TZ(2)
        COMMON /CCC6/ER(2),EQ(2),EZ(2),UV(2),VV(2),WV(2)
        COMMON /CCC7/CA(2),CB(2),CD(2),TD(2),TE(2),TF(2)
        COMMON /CCC8/EG(2),EH(2),EI(2),UH(2),VH(2),WH(2)
        COMMON /CCC9/CRQ,CQZ,CZR,TRQ,TQZ,TZR,US,VS,WS
        SN=SIN(CE)
        CS=COS(CE)
        S1=R1*SN
        R2=SQRT(R1*R1+D*D+2*D*S1)
	  IF(R2.EQ.0)THEN
	  S1=0
	  C1=-1
	  S2=-1
	  C2=0
	  ELSE
        S1=(S1+D)/R2
        C1=CS/R2
        S2=D*C1
        C1=R1*C1
        C2=(R1+D*SN)/R2
	  ENDIF
        F1=S2*S2
        F2=C2*C2
        F3=S2*C2
        F4=S2*S1
        F5=C2*C1
        F6=S2*C1
        F7=C2*S1
        IF(V.EQ.1) THEN
        P1=CR(1)+CR(2)*F2+CQ(2)*F1
        P2=CQ(1)+CR(2)*F1+CQ(2)*F2
        P3=CZ(1)+CZ(2)
        P4=(CR(2)-CQ(2))*F3
        P5=TZ(2)*S2
        P6=TZ(1)+TZ(2)*C2
        P7=UV(1)+UV(2)*C2
        P8=UV(2)*S2
        P9=WV(1)+WV(2)
        ENDIF
        IF(H.EQ.1) THEN
        PM=2*TD(2)*F3*S1
        Q1=CA(1)*CS+(CA(2)*F2+CB(2)*F1)*C1-PM
        Q2=CB(1)*CS+(CA(2)*F1+CB(2)*F2)*C1+PM
        Q3=CD(1)*CS+CD(2)*C1
        Q4=TD(1)*SN+(CA(2)-CB(2))*F3*C1+TD(2)*(F2-F1)*S1
        Q5=TE(1)*SN+TE(2)*F7+TF(2)*F6
        Q6=TF(1)*CS+TF(2)*F5-TE(2)*F4
        Q7=UH(1)*CS+UH(2)*F5-VH(2)*F4
        Q8=VH(1)*SN+UH(2)*F6+VH(2)*F7
        Q9=WH(1)*CS+WH(2)*C1
        ELSE
        DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/9*0/
        ENDIF
        CRQ=P1+F*Q1
        CQZ=P2+F*Q2
        CZR=P3+F*Q3
        TRQ=P4+F*Q4
        TQZ=P5+F*Q5
        TZR=P6+F*Q6
         US=P7+F*Q7
         VS=P8+F*Q8
         WS=P9+F*Q9
        RETURN
        END
********************************************
        !贝塞尔积分函数（）
	  UBROUTINE BESSEL(NT,Q,BL)
        DIMENSION BA(14),BB(14),BC(14)
        DATA BA/1.0,-2.2499997,1.26562079,-.3163866,.0444479,-.0039444,
     *  .00021,.5,-.56249985,.21093573,-.03954289,.00443319,
     *  -.00031761,.00001109/
        DATA BB/.79788456,-.00000077,-.0055274,-.00009512,.00137237,
     *  -.00072805,.00014476,.79788456,.00000156,.01659667,.00017105,
     *  -.00249511,.00113653,-.00020033/
        DATA BC/.78539816,.04166397,.00003954,-.00262573,.00054125,
     *  .00029333,-.00013558,.78539816,-.12499612,-.0000565,.00637879,
     *  -.00074348,-.00079824,.00029166/
        IF(Q.GE.3) GOTO 110
        BX=Q*Q/9
        BL=0
        DO 100 K=NT*7+7,NT*7+1,-1
100     BL=BL*BX+BA(K)
        IF(NT.EQ.1) BL=BL*Q
        RETURN
110     BX=3/Q
        BI=0
        BJ=0
        DO 120 K=NT*7+7,NT*7+1,-1
        BI=BI*BX+BB(K)
120     BJ=BJ*BX+BC(K)
        BL=COS(Q-BJ)
        IF(NT.EQ.1) BL=SIN(Q-BJ)
        BL=BL*BI/SQRT(Q)
        RETURN
        END
********************************************
        SUBROUTINE FFUNC(AQ,BQ,R,ZQ,MP,TR)
        TR=1
        O1=1
        O2=1
        DO 100 I=1,MP
        O1=O1*(AQ+I-1)/(R+I-1)
        O2=O2*(BQ+I-1)/I
100     TR=TR+O1*O2*ZQ**I
        RETURN
        END
*********************************************
        SUBROUTINE PRTRES(PO)
        COMMON /CCC9/CR,CQ,CZ,TR,TQ,TZ,US,VS,WS
        COMMON /CCC10/C1,C2,C3,E1,E2,E3,ER,EQ,EZ
        E1=TR*TR
        E2=TQ*TQ
        E3=TZ*TZ
        C1=CQ+CZ
        C2=CQ*CZ
        F1=CR+C1
        F2=CR*C1+C2-E1-E2-E3
        F3=CR*C2+2*TR*TQ*TZ-CR*E2-CQ*E3-CZ*E1
        E1=F1/3
        C1=F2-F1*E1
        C2=F2*E1-2*E1**3-F3
        E2=SQRT(-C1/3)
        Y=-0.5*C2/E2**3
        C1=ACOS(Y)
        C2=C1/3
        E3=COS(C2)
        C3=SQRT(3.)*SIN(C2)
        C1=E1+2*E2*E3
        C2=E1-E2*(E3-C3)
        C3=E1-E2*(E3+C3)
        E1=C1-PO*(C2+C3)
        E2=C2-PO*(C3+C1)
        E3=C3-PO*(C1+C2)
        ER=CR-PO*(CQ+CZ)
        EQ=CQ-PO*(CZ+CR)
        EZ=CZ-PO*(CR+CQ)
        RETURN
        END
******************************************
        SUBROUTINE MMPP(RD,MP)
        IF(RD.EQ.0)    MP=1
        IF(RD.GT.0     .AND.RD.LE.0.25)   MP=6
        IF(RD.GT.0.25  .AND.RD.LE.0.50)   MP=8
        IF(RD.GT.0.50  .AND.RD.LE.0.75)   MP=14
        IF(RD.GT.0.75  .AND.RD.LE.0.901)  MP=36
        IF(RD.GT.0.901 .AND.RD.LE.0.951)  MP=70
        IF(RD.GT.0.951 .AND.RD.LE.0.991)  MP=344
        IF(RD.GT.0.991 .AND.RD.LE.0.9951) MP=684
        IF(RD.GT.0.9951.AND.RD.LT.1.0)    MP=1000
        IF(RD.GT.1.0   .AND.RD.LE.1.009)  MP=1000
        IF(RD.GT.1.009 .AND.RD.LE.1.049)  MP=346
        IF(RD.GT.1.049 .AND.RD.LE.1.099)  MP=72
        IF(RD.GT.1.099 .AND.RD.LE.1.249)  MP=38
        IF(RD.GT.1.249 .AND.RD.LT.1.5)    MP=18
        IF(RD.GE.1.5   .AND.RD.LT.2.0)    MP=10
        IF(RD.GE.2.0   .AND.RD.LT.2.5)    MP=8
        IF(RD.GE.2.5   .AND.RD.LE.3.0)    MP=6
        IF(RD.GT.3.0)  MP=4
        RETURN
        END
**************************************************
        SUBROUTINE SURFA(M,PO,RR,YW,YI,H,V,SSU,SSV)
        DIMENSION SSU(5),SSV(8),RR(2),YW(2,4),YI(2,3)
        R6=1-2*PO
        R7=1-R6
        RD=RR(M)
        R8=RD*RD
        R9=1+PO
        IF(RD.EQ.1) GOTO 35
        CALL MMPP(RD,MP)
        IF(RD.NE.0) RA=1/R8
        IF(RD-1.) 30,35,40
30      CALL FFUNC(.5,-.5,1.,R8,MP,V3)
        V1=1
        V4=.5*RD
        GOTO 50
40      CALL FFUNC(.5,.5,2.,RA,MP,TR)
        V1=0
        V3=.5*TR/RD
        V4=.5/RD
        GOTO 50
35      V1=.5
        V3=.636619772
        V4=.5
50      Y1=V1-YW(M,1)
        Y3=V3-YW(M,3)
        Y4=V4-YW(M,4)
        IF(V.EQ.1) THEN
        V2=.5
        IF(RD.GT.1.) V2=.5/R8
        Y2=V2-YW(M,2)
        SSU(1)=-Y1+R6*Y2
        SSU(2)=-Y1*R7-R6*Y2
        SSU(3)=-Y1
        SSU(4)=-.5*R9*R6*Y4
        SSU(5)=R9*(R9+R6-1)*Y3
        ENDIF
        IF(H.EQ.1) THEN
        IF(RD-1.) 60,65,70
60      IF(RD.EQ.0) THEN
        H1=0
        H2=0
        ELSE
        CALL FFUNC(1.5,.5,2.,R8,MP,TR)
        H1=.5*RD*TR
        CALL FFUNC(1.5,.5,3.,R8,MP,TR)
        H2=.125*RD*TR
        ENDIF
        GOTO 80
70      CALL FFUNC(1.5,.5,2.,RA,MP,TR)
        AA=.5*RA
        H1=AA*TR
        CALL FFUNC(1.5,-.5,2.,RA,MP,TR)
        H2=AA*TR
        GOTO 80
65      H2=.21220659
        H1=0
80      H3=H2*RD
        X1=H1-YI(M,1)
        X2=H2-YI(M,2)
        X3=H3-YI(M,3)
        SSV(1)=2*X1+R7*X2
        SSV(2)=R7*(X1-X2)
        SSV(3)=R7*X2-X1
        SSV(4)=-Y1
        SSV(5)= Y1
        SSV(6)=-.25*R7*R9*X3-.5*(R6+R9)*R9*Y3
        SSV(7)=-.25*R7*R9*X3+.5*(R6+R9)*R9*Y3
        SSV(8)=-.5*R6*R9*Y4
        ENDIF
        RETURN
        END
************************************************
        SUBROUTINE SOEQ(L,N,A,Q)
        COMMON /CCC3/F1(24),F2(24),F3(24),F4(24)
        COMMON /CCC4/FA(24),FB(24),FC(24),FD(24),FE(24),FF(24)
        DOUBLE PRECISION A(144,12),Q(144),EZ,SS,T,R
        M=L*N
        DO 150 K=0,N-1
        L1=L*K
        M1=L1+L*3/2
        IF(K.EQ.N-1) M1=L1+L
        DO 100 I=L1+1,L1+L/2
        DO 100 J=1,L
        A(I,J)=A(I,J+L)
100     A(I,J+L)=0
        DO 150 I=L1+1,L1+L
        II=I-L1
        IP=0
        EZ=0
        DO 110 J=I,M1
        IF(ABS(A(J,II)).GT.ABS(EZ)) THEN
        EZ=A(J,II)
        IP=J
        ENDIF
110     CONTINUE
        IF(IP.EQ.I) GOTO 125
        DO 120 J=II,2*L
        SS=A(I,J)
        A(I,J)=A(IP,J)
120     A(IP,J)=SS
        T=Q(I)
        Q(I)=Q(IP)
        Q(IP)=T
125     R=1/A(I,II)
        DO 130 J=II+1,2*L
130     A(I,J)=R*A(I,J)
        Q(I)=R*Q(I)
        IF(I.EQ.M) GOTO 160
        DO 150 J=I+1,M1
        DO 140 JJ=II+1,2*L
140     A(J,JJ)=A(J,JJ)-A(J,II)*A(I,JJ)
150     Q(J)=Q(J)-A(J,II)*Q(I)
160     L3=L*(N-1)-1
        DO 170 I=M-1,M-L+1,-1
        DO 170 J=I,M-1
170     Q(I)=Q(I)-A(I,J-L3)*Q(J+1)
        DO 180 K=N-2,0,-1
        L2=L*K
        DO 180 I=L2+L,L2+1,-1
        II=I-L2
        DO 180 J=II,2*L-1
180     Q(I)=Q(I)-A(I,J+1)*Q(J+1+L2)
        IF(L.EQ.4) THEN
        DO 190 I=1,N
        J=4*I
        F3(I)=Q(J-3)
        F1(I)=Q(J-2)
        F4(I)=Q(J-1)
190     F2(I)=Q(J)
        ENDIF
        IF(L.EQ.6) THEN
        DO 200 I=1,N
        J=6*I
        FC(I)=Q(J-5)
        FA(I)=Q(J-4)
        FF(I)=Q(J-3)
        FD(I)=Q(J-2)
        FB(I)=Q(J-1)
200     FE(I)=Q(J)
        ENDIF
        RETURN
        END
************************ THE END *********************