开源Johnson-Cook损伤vumat子程序
2026年5月19日 21:39浏览：6
Johnson-CooK （简称 JC）模型主要用于解决金属材料在强冲击、高应变率、剧烈温度变化下的复杂响应问题。在国防穿甲爆破、航空航天器外壳受撞击、汽车高速碰撞以及工业上的金属切削加工等极端工况下，金属材料在极短时间内会发生巨大的变形，并且伴随着由于剧烈摩擦和变形产生的局部高温。传统的弹塑性模型无法准确模拟这种“又快、又热、变形又大”的极端物理过程，而 JC 模型正是为了破解这些高能耗、高破坏性的力学难题而诞生的。
该模型的核心思想是将复杂的金属材料行为进行“解耦”，认为材料的强度主要受到三个独立因素的叠加影响：应变硬化、应变率（变形速度）强化和热软化。简单来说，它认为金属材料在变形时有三个特点：一是随着变形量增大材料会越变越硬；二是变形发生得越快材料也会变得越硬；三是当变形产生的热量让材料温度升高时，材料就会变软。同时，模型还引入了热功转换机制，将材料变形产生的绝热塑性功直接转化为热量，并配合损伤退化和单元删除机制，从而能够逼真地模拟出材料从开始变形、变硬、变软，直到最终断裂撕裂的全过程。
它之所以成为高应变率仿真领域的“长青树”，主要原因有三点。首先是参数物理意义明确且极易获取，相比其他复杂的力学模型，JC 模型的参数可以通过标准的高速拉伸或霍普金森压杆（SHPB）试验轻松测得，工程实用性极高。其次是计算效率与数值稳定性极佳，它的数学形式简洁高效，非常适合显式动力学子程序（如 VUMAT）进行大规模并行计算，不易发生数值发散。最后是完美闭环了“力-热-损伤”的耦合，它不仅能算应力，还能同步算出温度升高以及材料的受损程度，在模拟金属穿透、飞溅、切屑形成等断裂失效行为时，具有无与伦比的仿真精度和视觉逼真度。
这里分享一个经典的vumat子程序，方便大家学习Johnson-Cook的相关理论模型：
原始链接：https://github.com/mauroarcidiacono/Abaqus-VUMAT-Johnson-Cook/tree/main
代码由Arcidiacono, Mauro F. and Rahimi, Salaheddin等人开发
! ########################################################################! 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
作者比较了使用代码计算和abaqus内置的Johnson-Cook模型计算响应的比较（与abaqus内置模型保持一致的精度）：
对Johnson-Cook建模感兴趣的可以下载了解，或者在这个代码基础上进行二次开发。相关代码和案例也上传的知识星球，也可以加入共同探讨交流这个Johnson-Cook模型代码：
技术邻APP
工程师必备
工程师必备
- 项目客服
- 培训客服
- 平台客服
TOP
