开源Johnson-Cook损伤vumat子程序

! ########################################################################! User subroutine to model the Johnson-Cook plasticity, damage and the ! Taylor-Quinney conversion of mechanical work into heat during plastic! deformation.!! Abaqus version: Abaqus 2022! Intel Fortran Compiler 2021.11! Visual Studio 2019!! Author: Mauro Francisco Arcidiacono! ########################################################################! ########################################################################! ! State Variable (SV) Definitions! SV1: initiation flag. If 0, first step, otherwise, 1.! SV2: effective plastic strain.! SV3: temperature.! SV4: Von Mises yield stress.! SV5: plastic strain increment.! SV6: parameter D of the Johnson-Cook damage model.! SV7: damage evolution parameter.! SV8: total number of iterations.! SV9: status of the element. If 1, the element is active, otherwise, the! element was deleted.! SV10 to SV16: stress tensor in the previous step.!! Note: this code was made to run the job using double precision due to! variable declaration inside the subroutines (real*8).!! ########################################################################

      subroutine vumat(! Read only (unmodifiable)variables -     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, jInfoArray,     2  stepTime, totalTime, dtArray, cmname, coordMp, charLength,     3  props, density, strainInc, relSpinInc,     4  tempOld, stretchOld, defgradOld, fieldOld,     5  stressOld, stateOld, enerInternOld, enerInelasOld,     6  tempNew, stretchNew, defgradNew, fieldNew,! Write only (modifiable) variables -     7  stressNew, stateNew, enerInternNew, enerInelasNew )!      include 'vaba_param.inc'      parameter (i_info_AnnealFlag = 1,      *     i_info_Intpt    = 2, ! Integration station number     *     i_info_layer  = 3, ! Layer number     *     i_info_kspt   = 4, ! Section point number in current layer     *     i_info_effModDefn = 5, ! =1 if Bulk/ShearMod need to be defined     *     i_info_ElemNumStartLoc   = 6) ! Start loc of user element number!      dimension props(nprops), density(nblock), coordMp(nblock,*),     1  charLength(nblock), dtArray(2*(nblock)+1), strainInc(nblock,ndir+nshr),     2  relSpinInc(nblock,nshr), tempOld(nblock),      3  stretchOld(nblock,ndir+nshr),     4  defgradOld(nblock,ndir+nshr+nshr),     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),     6  stateOld(nblock,nstatev), enerInternOld(nblock),     7  enerInelasOld(nblock), tempNew(nblock),     8  stretchNew(nblock,ndir+nshr),     8  defgradNew(nblock,ndir+nshr+nshr),     9  fieldNew(nblock,nfieldv),     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),     2  enerInternNew(nblock), enerInelasNew(nblock), jInfoArray(*)!      character*80 cmname      integer num_iter, max_num_iter, i, j      real*8 E, nu, A, B, n, m, Tm, Tr, C, epsilon_dot_zero, D1, D2,     1 D3, D4, D5, beta, Cp, D, pl_disp_failure, dmg_evol,     2 convergence_tolerance_factor, tolerance, mu, lambda,      3 eps, Temp, sigmaY, equiv_stress, pl_strain_inc, trace_strainInc,      4 eps_iter, equiv_stress_jc, f, pl_strain_inc_min,      5 pl_strain_inc_max, dWork, dPwork, rho, equiv_strain_fracture,     6 C_mat(ndir+nshr, ndir+nshr), stress_old(ndir+nshr),     7 trial_stress(ndir+nshr), dev_stress(ndir+nshr),      8 pl_strain_dir(ndir+nshr), C_PSD(ndir+nshr),      9 corrected_stress_iter(ndir+nshr)!      pointer (ptrjElemNum, jElemNum)      dimension jElemNum(nblock)!      lAnneal = jInfoArray(i_info_AnnealFlag)       iLayer = jInfoArray(i_info_layer)      kspt   = jInfoArray(i_info_kspt)      intPt  = jInfoArray(i_info_Intpt)      iUpdateEffMod = jInfoArray(i_info_effModDefn)      iElemNumStartLoc = jInfoArray(i_info_ElemNumStartLoc)      ptrjElemNum = loc(jInfoArray(iElemNumStartLoc))!
! ############################ Input parameters ############################
      ! Elastic properties      E = props(1)                          ! Young's Modulus      nu = props(2)                         ! Poisson's Ratio
      ! Johnson-Cook (JC) plasticity parameters      A = props(3)                          ! Parameter A      B = props(4)                          ! Parameter B      n = props(5)                          ! Parameter n       m = props(6)                          ! Parameter m      Tm = props(7)                         ! Melting temperature      Tr = props(8)                         ! Reference/Room temperature      C = props(9)                          ! Parameter C      epsilon_dot_zero = props(10)          ! Reference strain rate
      ! Johnson-Cook damage parameters      D1 = props(11)      D2 = props(12)      D3 = props(13)      D4 = props(14)      D5 = props(15)
      ! Plastic displacement at failure      pl_disp_failure = props(16)
      ! Taylor-Quinney (TQ) parameters      beta = props(17)                      ! Taylor-Quinney coefficient      Cp = props(18)                        ! Heat capacity      rho = props(19)                       ! Density
      ! Convergence tolerance for the increment of effective strain calculation      convergence_tolerance_factor = props(20)      tolerance = E*convergence_tolerance_factor
      ! Maximum number of iterations      max_num_iter = props(21)

! ############################ Elasticity Matrix ############################
      ! Lame parameters      ! Lame first parameter      lambda = (E*nu)/((1.d0 + nu)*(1.d0 - 2.0d0*nu))        ! Lame second parameter        mu = E/(2.0d0*(1.d0 + nu))
      ! Elasticity matrix (C_mat)      ! Linear elastic, homogeneous and isotropic material       ! Initialize the elasticity matrix to zero      C_mat = 0.d0
      do i = 1, ndir          do j = 1, ndir              if (i == j) then                  C_mat(i, j) = lambda + 2.d0*mu              else                  C_mat(i, j) = lambda              end if          end do       end do       do i = ndir + 1, ndir + nshr          C_mat(i, i) = 2.d0*mu      end do 

! ###########################################################################################################      ! ####################### Loop through each integration point to perform computations #######################! ###########################################################################################################  
      ! Global loop starting point -> loops through each integration point of the model      ! The index of the integration points is i      do i = 1, nblock
          ! ############################## Initial state ###############################             if (stateOld(i, 1) == 0.d0) then
              ! Compute Hooke's Law as a function of the strain increment              ! Trace calculation (EPSILON_kk)              trace_strainInc = sum(strainInc(i, 1:ndir))
              ! Calculate the direct components              stressNew(i, 1:ndir) = stressOld(i, 1:ndir) +       1        lambda*trace_strainInc + 2.d0*mu*strainInc(i, 1:ndir)
              ! Calculate the shear components              stressNew(i, ndir+1:ndir+nshr) = stressOld(i, ndir+1:ndir+nshr) +     1        2.d0*mu*strainInc(i, ndir+1:ndir+nshr)  
              stateNew(i, 1) = 1.d0           ! Initiation flag              stateNew(i, 2) = 0.d0           ! Equivalent plastic strain              stateNew(i, 3) = Tr             ! Initial temperature              stateNew(i, 4) = 0.d0           ! Yield stress              stateNew(i, 5) = 0.d0           ! Plastic strain increment              stateNew(i, 6) = 0.d0           ! Parameter D              stateNew(i, 7) = 0.d0           ! Damage evolution parameter              stateNew(i, 8) = 1              ! Total number of iterations              stateNew(i, 9) = 1              ! Element status
              ! Store the stress tensor in SDVs for damage evolution computation              do j = 1, ndir+nshr                  stateNew(i, j + 9) = stressNew(i, j)              end do
          ! ############################ From step 2 onward ############################              else
              eps = stateOld(i, 2)                                      ! Previous effective plastic strain              Temp = stateOld(i, 3)                                     ! Temperature              sigmaY = stateOld(i, 4)                                   ! Yield stress              D = stateOld(i, 6)                                        ! Parameter D              dmg_evol = stateOld(i, 7)                                 ! Damage evolution parameter              trial_stress = 0.d0                                       ! Trial total stress with increment
              ! stress_old stores the stress tensor in the previous step              if (D < 1.d0) then                  stress_old(1:ndir+nshr) = stressOld(i, 1:ndir+nshr)                     else                  stress_old(1:ndir+nshr) = stateOld(i, 10:ndir+nshr+10)                  end if
              ! Compute Hooke's Law as a function of the strain increment              ! Calculates the stress trial increment tensor assuming a pure elastic response              call elastic_stress(strainInc(i, 1:ndir+nshr), trial_stress,       1                            stress_old(1:ndir+nshr), lambda, mu, ndir, nshr) 
              ! Start the calculations to obtain the direction and magnitude of the              ! plastic strain increment              ! Direction = df/dSigma              ! Magnitude = dLambda (plastic multiplier)              ! dEpsilon^p = dLambda*df/dSigma = dp*3/2*DevStress/EquivalentStress              call equivalent_stress(equiv_stress, trial_stress,      1                               dev_stress, ndir, nshr)
              ! Plastic strain direction              if (equiv_stress /= 0) then                  pl_strain_dir = (3.d0*dev_stress)/(2.d0*equiv_stress)              else                  pl_strain_dir = 0              end if
              ! Elasticity matrix C times the plastic strain direction to calculate              ! later the plastic corrector              C_PSD = matmul(C_mat, pl_strain_dir)
              ! Initial value of the strain increment              pl_strain_inc = 0.d0                                      ! Initial plastic strain increment              pl_strain_inc_min = 0.d0                                  ! Minimum plastic strain increment              pl_strain_inc_max = equiv_stress/(2.d0*mu)                ! Maximum plastic strain increment

              !##################################################################################################              ! ######## Iteration to obtain the initial plastic strain increment (pl_strain_inc) starts ########              ! Return mapping algorithm
              num_iter = 0    ! number of iterations
              do while (num_iter <= max_num_iter) 
                  ! Iteration control                  num_iter = num_iter + 1                  if (num_iter == max_num_iter) then                      print*, 'ERROR - too many iterations | iter = ', num_iter                      call XPLB_EXIT                   end if
                  ! Iteration effective plastic strain = effective plastic strain + increment                  eps_iter = eps + pl_strain_inc                    ! Effective plastic strain increment rate                    eps_rate = pl_strain_inc/dtArray(1)
                  ! Compute the stress using the plastic corrector for this iteration                   ! (corrected_stress_iter)                  corrected_stress_iter = trial_stress - pl_strain_inc*C_PSD
                  ! Calculate the equivalent stress with the corrected stress in this                  ! iteration                  call equivalent_stress(equiv_stress, corrected_stress_iter,      1                                   dev_stress, ndir, nshr)
                  ! Calculate the Johnson-Cook equivalent stress                  call johnson_cook_plasticity(eps_iter, A, B, n, m,      1                                         Tm, Tr, Temp, C,     2                                         epsilon_dot_zero,      3                                         eps_rate, equiv_stress_jc)
                  ! If the calculated JC is less than the previous one, take the previous one                  ! as equivalent yield stress.                   if (equiv_stress_jc < sigmaY) then                      equiv_stress_jc = sigmaY                  end if
                  ! Compare the calculated JC stress with the equivalent stress resulting from                  ! the corrected stress state. These two values should be as close as possible                  ! within the specified tolerance to obtain the plastic strain increment for a                  ! JC material                  f = equiv_stress - equiv_stress_jc
                  if (abs(f) < tolerance) then                      exit                  end if
                  ! If there is no plastic strain increment and the JC yield stress is bigger                  ! than the equivalent stress, the stress state is elastic                  if ((pl_strain_inc == 0.d0) .and. (f < 0.d0)) then                      exit                  end if
                  ! Update the min and max plastic strain increment limits                  if ((f >= 0.d0) .and. (pl_strain_inc >= pl_strain_inc_min)) then                      pl_strain_inc_min = pl_strain_inc                    end if    
                  if ((f < 0.d0) .and. (pl_strain_inc < pl_strain_inc_max)) then                      pl_strain_inc_max = pl_strain_inc                    end if
                  ! Update the plastic strain increment                  pl_strain_inc = 0.5d0 * (pl_strain_inc_max + pl_strain_inc_min)
                  ! Update the pl_strain_inc to ensure continuity across increments and to                   ! improve the initial guess                  if (num_iter == 1) then                      pl_strain_inc = stateOld(i, 5)                  end if
              end do
              ! Save the newly calculated stress              stressNew(i, 1:ndir+nshr) = corrected_stress_iter(1:ndir+nshr)
              ! Store the stress tensor in SDVs for damage evolution computation              do j = 1, ndir+nshr                  stateNew(i, j + 9) = stressNew(i, j)              end do
              ! Calculate the work increment              dWork = dot_product(0.5d0 * (corrected_stress_iter(1:ndir+nshr) + stress_old(1:ndir+nshr)),      1                            strainInc(i, 1:ndir+nshr))
              ! Calculate the plastic work increment              dPwork = 0.5d0 * pl_strain_inc * equiv_stress
              ! Calculate the internal energy per unit mass              enerInternNew(i) = enerInternOld(i) + dWork/rho
              ! Calculate the dissipated inelastic energy per unit mass              enerInelasNew(i) = enerInelasOld(i) + dPwork/rho
              ! Evaluate the equivalent strain at fracture              call johnson_cook_damage(equiv_strain_fracture, equiv_stress, Tm, Tr, D1, D2,     1                                 D3, D4, D5, Temp, epsilon_dot_zero, eps_rate,      2                                 corrected_stress_iter, ndir, nshr)
              ! Update the parameter D of the Johnson-Cook damage model              if (equiv_strain_fracture /= 0.d0) then                  D = D + abs(pl_strain_inc / equiv_strain_fracture)              end if
              ! Fracture is allowed to occur when D = 1.0              if (D >= 1.d0) then                  D = 1.d0                  dmg_evol = dmg_evol +      1            abs(pl_strain_inc * charLength(i) / pl_disp_failure)                  if (dmg_evol >= 1.d0) then                      dmg_evol = 1.d0                  end if                  stressNew(i, 1:ndir+nshr) = (1.d0 - dmg_evol)*stressNew(i, 1:ndir+nshr)              end if
              ! Update the state variables              stateNew(i, 1) = 1.d0                                                           ! Initiation flag              stateNew(i, 2) = eps_iter                                                       ! Equivalent plastic strain              stateNew(i, 3) = Temp + beta*dPwork/rho/Cp                                      ! Temperature              stateNew(i, 4) = equiv_stress_jc                                                ! Yield stress              stateNew(i, 5) = pl_strain_inc                                                  ! Plastic strain increment in the step              stateNew(i, 6) = D                                                              ! Parameter D              stateNew(i, 7) = dmg_evol                                                       ! Damage evolution parameter              stateNew(i, 8) = stateOld(i, 8) + num_iter                                      ! Total number of iterations
              if (dmg_evol >= 1.d0) then                  stateNew(i, 9) = 0    ! Delete element              else                  stateNew(i, 9) = 1    ! Active element              end if
          end if
      end do
      return      end
! Additional subroutines      include 'utils.for'       subroutine elastic_stress (strainInc, stressNew, stressOld,     1 lambda, mu, ndir, nshr)      ! This subroutine calculates the elastic stress tensor in Voigt      ! notation using a strain tensor. The calculation is done according       ! to Hooke's law.      ! SIGMA_ij = lambda*EPSILON_kk*KRONECKERDELTA_ij + 2*mu*EPSILON_ij
          integer ndir, nshr          real*8 strainInc(ndir+nshr), stressNew(ndir+nshr),       1    stressOld(ndir+nshr), lambda, mu, trace_strainInc 
          ! Trace calculation (EPSILON_kk)          trace_strainInc = sum(strainInc(1:ndir))
          ! Calculate the direct components          stressNew(1:ndir) = stressOld(1:ndir) +       1    lambda*trace_strainInc + 2.d0*mu*strainInc(1:ndir)
          ! Calculate the shear components          stressNew(ndir+1:ndir+nshr) = stressOld(ndir+1:ndir+nshr) +     1    2.d0*mu*strainInc(ndir+1:ndir+nshr)
      return      end

      subroutine equivalent_stress (equiv_stress, stress_tensor,      1 dev_stress, ndir, nshr)      ! This subroutine calculates the equivalent Von Mises stress             ! from a Voigt notation stress tensor.
          integer ndir, nshr          real*8 stress_tensor(ndir+nshr), hyd,      1    dev_stress(ndir+nshr), equiv_stress
          ! Calculate the hydrostatic component of the stress tensor          hyd = sum(stress_tensor(1:ndir))/3.d0          ! Calculate the deviatoric tensor          dev_stress(1:ndir) = stress_tensor(1:ndir) - hyd          dev_stress(ndir+1:ndir+nshr) = stress_tensor(ndir+1:ndir+nshr)
          ! Compute the equivalent Von Mises stress          ! stressVM = sqrt(3/2 * dev_stress : dev_stress)          equiv_stress = sqrt(3.d0/2.d0 *(dev_stress(1)**2.d0 +      1    dev_stress(2)**2.d0 + dev_stress(3)**2.d0 +      2    2.d0*dev_stress(4)**2.d0 + 2.d0*dev_stress(5)**2.d0 +          3    2.d0*dev_stress(6)**2.d0))
      return      end

      subroutine johnson_cook_plasticity(eps_iter, A, B, n, m, Tm, Tr,      1 T, C, epsilon_dot_zero, eps_rate, equiv_stress_jc)      ! This subroutine calculates the Von Mises stress using the       ! Johnson Cook equation. 
          real*8 eps_iter, A, B, n, m, Tm, Tr, T, C,       1    epsilon_dot_zero, eps_rate, homologous_Temp,      2    equiv_stress_jc  
          if (T < Tr) then              homologous_Temp = 0.d0          else if (T > Tm) then              homologous_Temp = 1.d0          else              homologous_Temp = (T - Tr)/(Tm - Tr)          end if
          if (eps_rate <= 1D-12) then              eps_rate = epsilon_dot_zero          end if
          if (eps_iter <= 1D-12) then              eps_iter = 0.d0          end if
          equiv_stress_jc = (A + B*eps_iter**n) *     1    (1 + C*log(eps_rate/epsilon_dot_zero)) *     2    (1 - homologous_Temp**m)
      return      end

      subroutine johnson_cook_damage(equiv_strain_frac, equiv_stress,     1 Tm, Tr, D1, D2, D3, D4, D5, T, epsilon_dot_zero, eps_rate,     2 stress_tensor, ndir, nshr)      ! This subroutine calculates the equivalent strain for fracture      ! for the Johnson-Cook damage model. 
          integer ndir, nshr          real*8 Tm, Tr, D1, D2, D3, D4, D5, T,       1    epsilon_dot_zero, eps_rate, homologous_Temp,      2    pressure_stress_ratio, equiv_strain_frac,     3    stress_tensor(ndir+nshr), hyd, equiv_stress
          if (T < Tr) then              homologous_Temp = 0.d0          else if (T > Tm) then              homologous_Temp = 1.d0          else              homologous_Temp = (T - Tr)/(Tm - Tr)          end if
          if (eps_rate <= 1D-12) then              eps_rate = epsilon_dot_zero          end if
          ! Calculate the hydrostatic component of the stress tensor          hyd = sum(stress_tensor(1:ndir))/3.d0
          ! Calculate the pressure stress ratio          if (equiv_stress == 0.d0) then              pressure_stress_ratio = 0.d0          else              pressure_stress_ratio = hyd/equiv_stress          end if
          equiv_strain_frac = (D1 +      1    D2*exp(-D3*pressure_stress_ratio))*     2    (1 + D4*log(eps_rate/epsilon_dot_zero))*     3    (1 + D5*homologous_Temp)
      return      end