ammarkh95 / abaqus_pdalac Goto Github PK

  SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
 1 RPL,DDSDDT,DRPLDE,DRPLDT,
 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
 

  INCLUDE 'ABA_PARAM.INC' 

  CHARACTER*80 CMNAME
  DIMENSION STRESS(NTENS),STATEV(NSTATV),
 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
 3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
 4 JSTEP(4),UPSTRAN(NTENS),CD(NTENS,NTENS),Dgrd(6),dmg(4),Xphi(6),n(6)

  PARAMETER(ZERO = 0.0D0, ONE=1.0D0, TWO=2.0D0, PI = 3.141592653589793D0)
  
  INTEGER fflags(2)
    REAL(8) :: a, theta, theta_constant, T0, T90, delta, V21, V31, V32
    REAL(8) :: Q11, Q12, Q13, Q22, Q23, Q33, Q44, Q55, Q66
    REAL(8) :: cos_theta, sin_theta, cos2, sin2, cos4, sin4
    REAL(8) :: sin_cos3, cos_sin3
    REAL(8) :: z_coord
    REAL(8) :: layer_thickness, layer1_max, layer2_max, layer3_max, layer4_max
  

  layer_thickness = 0.025D0
  layer1_max = layer_thickness
  layer2_max = 2.0D0 * layer_thickness
  layer3_max = 3.0D0 * layer_thickness
  layer4_max = 4.0D0 * layer_thickness

  E1=PROPS(1)
  E2=PROPS(2)
  E3=PROPS(3)
  V12=PROPS(4)
  V13=PROPS(5)
  V23=PROPS(6)   
  G12=PROPS(7)
  G13=PROPS(8)
  G23=PROPS(9)

  DO I=1,NTENS
   DO J=1,NTENS
   DDSDDE(I,J)=0.0D0       
   ENDDO
  ENDDO

  V21 = V12 * E2 / E1
  V31 = V13 * E3 / E1
  V32 = V23 * E3 / E2

  delta = (1.0D0 - V12 * V21 - V23 * V32 - V13 * V31 - 2.0D0 * V21 * V32 * V13) / (E1 * E2 * E3)

  Q11 = (1.0D0 - V23 * V32) / (E2 * E3 * delta)
  Q12 = (V21 + V23 * V31) / (E2 * E3 * delta)
  Q13 = (V31 + V21 * V32) / (E2 * E3 * delta)
  Q22 = (1.0D0 - V13 * V31) / (E1 * E3 * delta)
  Q23 = (V32 + V12 * V31) / (E1 * E3 * delta)
  Q33 = (1.0D0 - V12 * V21) / (E1 * E2 * delta)
  Q44 = G23
  Q55 = G13
  Q66 = G12

  a = 1.0D0  
  
  T0 = 0.0D0 * (PI / 180.0D0)
  T90 = 90.0D0 * (PI / 180.0D0)

  z_coord = COORDS(3)

  IF (z_coord < layer1_max) THEN
      theta_constant = (2.0D0 * (T90 - T0) / a) * DABS(z_coord - a/2.0D0) + T0
      theta = theta_constant

  ELSE IF (z_coord < layer2_max) THEN
      theta_constant = (2.0D0 * (T90 - T0) / a) * DABS(z_coord - a/2.0D0) + T0
      theta = -theta_constant

  ELSE IF (z_coord < layer3_max) THEN
      theta_constant = (2.0D0 * (T90 - T0) / a) * DABS(z_coord) + T0
      theta = theta_constant

  ELSE IF (z_coord < layer4_max) THEN
      theta_constant = (2.0D0 * (T90 - T0) / a) * DABS(z_coord) + T0
      theta = -theta_constant

  END IF
    
  cos_theta = COS(theta)
  sin_theta = SIN(theta)

  cos2 = cos_theta * cos_theta
  sin2 = sin_theta * sin_theta

  cos4 = cos2 * cos2
  sin4 = sin2 * sin2

  sin_cos3 = sin_theta * cos2 * cos_theta
  cos_sin3 = cos_theta * sin2 * sin_theta
                 
  DDSDDE(1, 1) = Q11 * cos4 + Q22 * sin4 + (2.0D0 * (sin2 * cos2) * (Q12 + 2.0D0 * Q66))
                 
  DDSDDE(1, 2) = ((Q11 + Q22 - 4.0D0 * Q66) * (sin2 * cos2)) + (Q12 * (sin4 + cos4))
                 
  DDSDDE(1, 3) = Q13 * cos2 + Q23 * sin2
                 
  DDSDDE(1, 6) = (Q11 - Q12 - 2.0D0 * Q66) * sin_cos3 - (Q22 - Q12 - 2.0D0 * Q66) * cos_sin3

  DDSDDE(2, 1) = DDSDDE(1, 2)
                 
  DDSDDE(2, 2) = Q11 * sin4 + (2.0D0 * (Q12 + 2.0D0 * Q66)) * (sin2 * cos2) + Q22 * cos4
                 
  DDSDDE(2, 3) = Q13 * sin2 + Q23 * cos2
                 
  DDSDDE(2, 6) = (Q11 - Q12 - 2.0D0 * Q66) * cos_sin3 - (Q22 - Q12 - 2.0D0 * Q66) * sin_cos3

  DDSDDE(3, 1) = DDSDDE(1, 3)

  DDSDDE(3, 2) = DDSDDE(2, 3)

  DDSDDE(3, 3) = Q33
                 
  DDSDDE(3, 6) = (Q13 - Q23) * cos_theta * sin_theta
                 
  DDSDDE(4, 4) = Q44 * cos2 + Q55 * sin2
                    
  DDSDDE(4, 5) = (Q55 - Q44) * cos_theta * sin_theta

  DDSDDE(5, 4) = DDSDDE(4, 5)
                 
  DDSDDE(5, 5) = Q55 * cos2 + Q44 * sin2

  DDSDDE(6, 1) = DDSDDE(1, 6)

  DDSDDE(6, 2) = DDSDDE(2, 6)

  DDSDDE(6, 3) = DDSDDE(3, 6)
                    
  DDSDDE(6, 6) = (Q11 - 2.0D0 * Q12 + Q22) * (sin2 * cos2) + Q66 * ((sin2 - cos2)*(sin2 - cos2))

   DO I=1,NTENS  
	  DO J=1, NTENS
    
 	   STRESS(I)=STRESS(I)+DDSDDE(I,J)*DSTRAN(J)

      END DO
  END DO	  

  RETURN
  END