[其他行业算法] 关于滑动轴承承载力的求解

PROGRAM GASJOURNAL
    DIMENSION H(61,61),P(61,61),Y(61),X(61)
    DATA B,R,C0,AN,PA,EDA,EPSON/60.0E-3,25.0E-3,5.0E-5,6.0E4,1.013E5,1.79E-5,0.8/
    OPEN(8,FILE='PRESSURE.DAT',STATUS='UNKNOWN')
    OPEN(9,FILE='HIGHT.DAT',STATUS='UNKNOWN')
    OPEN(11,FILE='LOAD.DAT',STATUS='UNKNOWN')
    PI=3.1415926
    N=61
    N1=N-1
    M=61
    M1=M-1
    DX=2.0*PI/N1
    DY=2.0/M1
    OMEGA=AN*2.0*PI/60
    U=OMEGA*R
    ALENDA=6.0*EDA*U*R/PA/C0**2
    ALFA=(R/B*DX/DY)**2
    SUM1=0.0
    S=D*B

    CALL SUBH (N,M,DX,EPSON,H)
    CALL SUBP (N,M,DX,EPSON,ALFA,ALENDA,H,P)
    CALL SUBM(N,M,P)
    ALENDA=6.0*EDA*U*R/PA/C0**2
    SUM1=ALENDA*S*SUM1
    DO I=1,N
        X(I)=(I-1)*DX
    ENDDO

    DO J=1,M
        Y(J)=(J-1.0)*DY-0.5
    ENDDO

    WRITE(8,40)Y(1),(Y(J),J=1,M)
    DO I=1,N
        WRITE(8,40)X(I)*180.0/PI,(P(I,J),J=1,M)
    ENDDO
    WRITE(9,40)Y(1),(Y(J),J=1,M)
        DO I=1,N
        WRITE(9,40)X(I)*180.0/PI,(H(I,J),J=1,M)
    END DO
!WRITE(8,40)AX,(H(I,J),J=1,M)
!WRITE(9,40)AX,(P(I,J),J=1,M)
    WRITE(11,*) '承载力=',SUM1
40  FORMAT(62(E12.6,1X))
    STOP
    END

SUBROUTINE SUBH(N,M,DX,EPSON,H)
    DIMENSION H(N,M)
        DO I=1,N
        SETA=(I-1.0)*DX
        DO J=1,M
        H(I,J)=1.0+EPSON*COS(SETA)
        ENDDO
        ENDDO
        RETURN
    END

    SUBROUTINE SUBP(N,M,DX,EPSON,ALFA,ALENDA,H,P)
        DIMENSION H(N,M),P(N,M)
        DO I=2,N-1
        DO J=2,M-1
        P(I,J)=1.1
        ENDDO   
        ENDDO
        DO I=1,N
        P(I,1)=1.0  
        P(I,M)=1.0
        ENDDO
        DO J=1,M
        P(1,J)=1.0
        P(N,J)=1.0
        ENDDO
        IK=0
10        C1=0.0
        ALOAD=0.0
        DO I=2,N-1
        I1=I-1
        I2=I+1
        DO J=2,M-1
        J1=J-1
        J2=J+1
        PD=P(I,J)
    A1=H(I,J)**3+(3*H(I,J)**2*(H(I2,J)-H(I1,J))/4)
    A2=H(I,J)**3-(3*H(I,J)**2*(H(I2,J)-H(I1,J))/4)
    A3=2*H(I,J)**3*(1+ALFA)
    A4=ALFA*(H(I,J)**3+(3*H(I,J)**2*(H(I,J2)-H(I,J1))/4))
    A5=ALFA*(H(I,J)**3-(3*H(I,J)**2*(H(I,J2)-H(I,J1))/4))
    A6=(ALENDA*DX/2)*(P(I,J)*(H(I2,J)-H(I1,J))+ H(I,J)*(P(I2,J)-P(I1,J)))
    P(I,J)=(A1*P(I2,J)**2+A2*P(I1,J)**2+A4*P(I,J2)**2+A5*P(I,J1)**2-A6)/A3
    !A1=(0.5*(H(I2,J)+H(I,J)))**3
        !A2=(0.5*(H(I,J)+H(I1,J)))**3
    !A3=ALFA*(0.5*(H(I,J2)+H(I,J)))**3
        !A4=ALFA*(0.5*(H(I,J)+H(I,J1)))**3
        !P(I,J)=(-DX*ALENDA*(P(I+1,J)*H(I+1,J)-P(I-1,J)*H(I-1,J))+A1*P(I2,J)**2+A2*P(I1,J)**2+A3*P(I,J2)**2+A4*P(I,J1)**2)/(A1+A2+A3+A4)
        IF(P(I,J).LT.0.0)P(I,J)=0.0
        P(I,J)=SQRT(P(I,J))
        P(I,J)=0.7*PD+0.3*P(I,J)
        IF(P(I,J).LT.0.0)P(I,J)=0.0
        C1=C1+ABS(P(I,J)-PD)
        ALOAD=ALOAD+P(I,J)
        ENDDO
        ENDDO
        IK=IK+1
        C1=C1/ALOAD
        WRITE(*,*)IK,C1,ALOAD
        IF(C1.GT.1.E-7)GOTO 10
        RETURN
    END

    SUBROUTINE SUBM(N,M,P)
    DIMENSION P(N,M)
    PX=0
    PY=0
    DO I=1,N
        AI=(I-1)*DX
        DO J=1,M
            PX=PX-(P(I,J)-1)*COS(AI)*DX*DY
            PY=PY+(P(I,J)-1)*SIN(AI)*DX*DY
        ENDDO
    ENDDO
    SUM1=SQRT(PX*PX+PY*PY)
    RETURN
    END