開源Johnson-Cook損傷vumat子程序
瀏覽：846 收藏：1
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
1
