mainstrs.f90

非线性有限元分析程序

        SUBROUTINE MainStrs( vstrs, vmain, itype )
!........
!       模块功能
!           计算主应力
!........
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION vstrs( 6 ), vmain( 3 ), vwork( 6 )
        vmain( 1 ) = 0.0D0
        vmain( 2 ) = 0.0D0
        vmain( 3 ) = 0.0D0
        IF( itype .EQ. 1 ) THEN
!.........单向应力状态
          vmain( 1 ) = vstrs( 1 )
        ELSE IF( itype .EQ. 2 ) THEN
!.........平面应力/应变状态
          cnst1 = ( vstrs( 1 ) + vstrs( 2 ) ) / 2.0D0
          cnst2 = ( vstrs( 1 ) - vstrs( 2 ) ) / 2.0D0
          cnst2 = DSQRT( cnst2 * cnst2 + vstrs( 3 ) * vstrs( 3 ) )
          vmain( 1 ) = cnst1 + cnst2
          vmain( 2 ) = cnst1 - cnst2
          vmain( 3 ) = vstrs( 4 )
        ELSE IF( itype .EQ. 3 ) THEN
!.........空间应力状态
          DO istrs = 1, 6
            vwork( istrs ) = vstrs( istrs )
          END DO
          cnst1 = ( vstrs( 1 ) + vstrs( 2 ) + vstrs( 3 ) ) / 3.0D0
          DO istrs = 1, 3
            vwork( istrs ) = vwork( istrs ) - cnst1
            vmain( istrs ) = cnst1
          END DO
          cnst2 = vwork( 1 ) * vwork( 2 ) + vwork( 2 ) * vwork( 3 ) +          &
                  vwork( 3 ) * vwork( 1 ) - vwork( 4 ) * vwork( 4 ) -          &
                  vwork( 5 ) * vwork( 5 ) - vwork( 6 ) * vwork( 6 )
          cnst3 = vwork( 1 ) * vwork( 2 ) * vwork( 3 ) +                       &
                  vwork( 4 ) * vwork( 5 ) * vwork( 6 ) * 2.0D0 -               &
                  vwork( 1 ) * vwork( 4 ) * vwork( 4 ) -                       &
                  vwork( 2 ) * vwork( 5 ) * vwork( 5 ) -                       &
                  vwork( 3 ) * vwork( 6 ) * vwork( 6 )
          IF( cnst2 .LT.-1.0D-8 ) THEN
            cnst2 = DSQRT(-cnst2 / 3.0D0 )
            angle = cnst3 / ( 2.0D0 * cnst2 * cnst2 * cnst2 )
            IF( angle .GT. 1.0D0 ) angle = 1.0D0
            IF( angle .LT.-1.0D0 ) angle =-1.0D0
            angle = ACOS( angle ) / 3.0D0
            cnst2 = 2.0D0 * cnst2
            vmain( 1 ) = vmain( 1 ) + cnst2 * DCOS( angle )
            vmain( 2 ) = vmain( 2 ) + cnst2 * DCOS( angle + 2.0943951D0 )
            vmain( 3 ) = vmain( 3 ) + cnst2 * DCOS( angle + 4.1887902D0 )
          END IF
        ELSE IF( itype .EQ. 4 ) THEN
!.........轴对称应力状态
          cnst1 = ( vstrs( 1 ) + vstrs( 3 ) ) / 2.0D0
          cnst2 = ( vstrs( 1 ) - vstrs( 3 ) ) / 2.0D0
          cnst2 = DSQRT( cnst2 * cnst2 + vstrs( 4 ) * vstrs( 4 ) )
          vmain( 1 ) = cnst1 + cnst2
          vmain( 2 ) = cnst1 - cnst2
          vmain( 3 ) = vstrs( 2 )
          IF( vmain( 2 ) .LT. vmain( 3 ) ) THEN
            vmain( 2 ) = vstrs( 2 )
            vmain( 3 ) = cnst1 - cnst2
          END IF
        END IF
        IF( itype .GE. 3 ) THEN
          IF( vmain( 1 ) .LT. vmain( 2 ) ) THEN
            angle = vmain( 1 )
            vmain( 1 ) = vmain( 2 )
            vmain( 2 ) = angle
          END IF
          IF( vmain( 1 ) .LT. vmain( 3 ) ) THEN
            angle = vmain( 1 )
            vmain( 1 ) = vmain( 3 )
            vmain( 3 ) = angle
          END IF
          IF( vmain( 2 ) .LT. vmain( 3 ) ) THEN
            angle = vmain( 2 )
            vmain( 2 ) = vmain( 3 )
            vmain( 3 ) = angle
          END IF
        END IF
        RETURN
        END