UMAT子程序详解

c ************************************************************************************** ************************************************************************************** **计算材料阵D **通过输入的Engineering Constants参数，计算材料弹性阵D。由于D为对称矩阵，只需计算9个值， **再通过对称赋值即可。 ************************************************************************************** C FILL THE 6X6 FULL STIFFNESS MATRIX DO I = 1, 6 DO J = 1, 6 CFULL(I,J)=ZERO END DO END DO ATEMP = ONE - TWO * XNULT * XNUTL - XNUTT ** TWO 1 - TWO * XNULT * XNUTL * XNUTT CFULL(1,1) = TENL * (ONE - XNUTT ** TWO) / ATEMP CFULL(2,2) = TENT * (ONE - XNULT * XNUTL) / ATEMP CFULL(3,3) = CFULL(2,2) CFULL(1,2) = TENT * (XNULT + XNULT * XNUTT) / ATEMP CFULL(1,3) = CFULL(1,2) CFULL(2,3) = TENT * (XNUTT + XNULT * XNUTL) / ATEMP CFULL(4,4) = SHRLT CFULL(5,5) = SHRLT CFULL(6,6) = SHRTT DO I = 2, 6 DO J = 1, I-1 CFULL(I,J) = CFULL(J,I) END DO END DO ************************************************************************************** ************************************************************************************** **计算失效应变值 **由于材料初始失效准则及损伤演化准则的要求，需要计算失效应变值。EPITL = σf,t L �C11 f,t f,c f ⁄ **EPICL = σf,c � C EPITT = σ � C EPICT = σ � C EPISLT = τ C 11 22 22 44 LT L T T ************************************************************************************** c calculate the failure strain by failure stress EPITL = SIGTL / cfull(1,1) !FAILURE STRAIN 1 DIRECTION IN TENSION EPICL = SIGCL / cfull(1,1) !FAILURE STRAIN 1 DIRECTION IN COMPRESSION EPITT = SIGTT / cfull(2,2) !TENSILE FAILURE STRAIN 2 DIRECTION EPICT = SIGCT / cfull(2,2) !COMPRESSIVE FAILURE STRAIN 2 DIRECTION EPISLT = SIGSLT/ SHRLT ! FAILURE SHEAR STRAIN ...ENGINEERING STRAIN ************************************************************************************** ************************************************************************************** **计算初始失效状态值 **调用函数CheckFailureIni()函数计算各个积分点初始失效状态值， 用于判断各个积分点是否进入 **失效状态,具体失效状态值存于STATEV中。
**** ***************************************************************************** **UMAT子程序接口，接口中的全局变量可以参考Abaqus User Subroutines Reference Manual **中，关于UMAT的介绍，对所有全局变量参数给出详细的介绍 ********************************************************************************* aC 4X4 JACOBIAN THEN CONDENSATION WITH REGULARIZATION c (energy due to viscous regularization is calculated) 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,KSTEP,KINC) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 CMNAME DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS), 2 DDSDDT(NTENS),DRPLDE(NTENS), 3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) ********************************************************************************** ********************************************************************************** **定义子程序局部变量，用于子程序内部计算使用。 **下面对各个局部变量给出英文解释 ********************************************************************************** DIMENSION STRANT(6),TSTRANT(4) DIMENSION CFULL(6,6),CDFULL(6,6) DIMENSION DDFDE(6), DDMDE(6), DCDDF(6,6), DCDDM(6,6) DIMENSION ATEMP1(6), ATEMP2(6), TDDSDDE(6,6) DIMENSION OLD_STRESS(6) DIMENSION DOLD_STRESS(6),D_STRESS(6) PARAMETER (ZERO = 0.D0,ONE = 1.D0,TWO = 2.D0, HALF = 0.5D0) C**************************** C STRANT..... STRAIN AT THE END OF THE INCREMENT C TSTRANT.....TEMPORARY ARRAY TO HOLD THE STRAIN FOR PLANE STRESS PROBLEM C CFULL.......FULL 6X6 ELASTICITY MATRIX C CDFULL......FULL 6X6 DAMAGED ELASTICITY MATRIX C DDFDE....... D DF/D E C DDMDE....... D DM/D E C DCDDF....... D C/ D DF THE DERIVATIVE OF THE FULL MATRIX OVER DF C DCDDM........D C/ D DM THE DERIVATIVE OF THE FULL MATRIX OVER DM C ATEMP1,ATEMP2...TEMPORARY ARRAY USED IN JACOBIAN CALCULATION C TDDSDDE.....UNCONDENSED JACOBIAN MATRIX FOR PLANE STRESS PROBLEM C OLD_STRESS...STRESS AT THE BEGINNING OF THE INCREMENT, SAVED FOR THE ENERGY C COMPUTATION C DOLD_STRESS...STRESS AT THE BEGINNING OF THE INCREMENT, C IF THERE'S NO VISCOUS REGULARIZATION C D_STRESS...STRESS IF THERE'S NO VISCOUS REGULARIZATION, THE ABOVE IS CALCULATED
C TO CALCULATE THE SCD, ENERGY CAUSED BY VISCOUS REGULARIZATION C STATEV(1) damage variable df C STATEV(2) damage variable dm C STATEV(3) regularized damage variable dfv C STATEV(4) regularizaed damage variable dmv C STATEV(5:10) TEMPORARY ARRAYS TO SAVE DOLD_STRESS C************************** ************************************************************************************** ************************************************************************************** **在INP文件中，调用*USER MATERIAL,将所输入材料参数赋予局部变量矩阵PROPS(14) **PROPS(1)--PROPS(6)为使用Engineering Constants给出的材料参数； **PROPS(7)—PROPS(14)为所定义材料失效参数、断裂能参数及粘性参数。 ************************************************************************************** C C GET THE MATERIAL PROPERTIES---ENGINEERING CONSTANTS C TENL = PROPS(1) !YOUNG'S MODULUS IN DIRECTION 1 (L) TENT = PROPS(2) !YOUNG'S MODULUS IN DIRECTION 2 (T) SHRLT = PROPS(3) !SHEAR MODULUS IN 12 PLANE SHRTT = PROPS(4) !SHEAR MODULUS IN 23 PLANE XNULT = PROPS(5) !POISON'S RATIO POI_12 XNUTT = PROPS(6) !POISON'S RATIO POI_23 μ μ XNUTL = XNULT / TENL * TENT !POI_21 由关系 TL = LT 得到 C C C GET THE FAILURE PROPERTIES
ET EL
SIGTL = PROPS(7) !FAILURE STRESS IN 1 DIRECTION IN TENSION SIGCL = PROPS(8) !FAILURE STRESS IN 1 DIRECTION IN COMPRESSION SIGTT = PROPS(9) !FAILURE STRESS IN 2 DIRECTION IN TENSION SIGCT = PROPS(10) !FAILURE STRESS IN 2 DIRECTION IN COMPRESSION SIGSLT = PROPS(11) !FAILURE STRESS IN SHEAR IN 1-2 PLANE GFMAT = PROPS(12) !FRACTURE ENERGY IN MATRIX GFFIB = PROPS(13) !FRACTURE ENERGY IN FIBER ETA = PROPS(14) ! VISCOSITY FOR REGULARIZATION ************************************************************************************** ************************************************************************************** **更新初始状态的应变分量 **STRAN为第K增量步初始应变，DSTRAN为上次迭代步结束时应变增量，则STRANT为当前迭代步初始 **应变量。在同一迭代步中，STRAN不变，DSTRAN应变增量不断变化，直到迭代收敛时，将 **STRANT=STRAN+DSTRAN存入第K+1步的STRAN，作为第K+1步的初始应变量。 ************************************************************************************** C C CALCULATE THE STRAIN AT THE END OF THE INBaidu NhomakorabeaREMENT C DO I = 1, NTENS STRANT(I) = STRAN(I) + DSTRAN(I) END DO