你帮我检查一下我的代码 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)
C-----------------------------------------------------------------------
C ABAQUS UMAT - 混凝土弹塑性本构模型（修复后版本）
C 修复内容：
C 1. 剪切模量EG初始化及传递问题
C 2. 拉伸软化区R_THETA_NR未赋值问题
C 3. EPSILON_V_YIELD更新逻辑
C 4. dETA_dLAMBDA默认初始化
C-----------------------------------------------------------------------
INCLUDE 'ABA_PARAM.INC'
CHARACTER*80 CMNAME
INTEGER :: NDI, NSHR, NTENS, NSTATV, NPROPS, NOEL, NPT, LAYER,
1 KSPT, KINC, JSTEP(4)
DOUBLE PRECISION :: 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 CELENT, PNEWDT, SSE, SPD, SCD, RPL, DRPLDT, DTEMP, TEMP
! --- 变量声明 ---
INTEGER :: I, J, N_SUBSTEPS, SUBSTEP, ITER
INTEGER :: N_ETA_PAIRS, I_ETA
DOUBLE PRECISION :: DT_FACTOR, TOTAL_STRAIN, MAX_PLASTIC_INCR
DOUBLE PRECISION :: STRESS_SUB(NTENS), EPLAS_SUB(NTENS)
DOUBLE PRECISION :: DSTRAN_SUB(NTENS), STRESS_TRIAL_SUB(NTENS)
DOUBLE PRECISION :: ERROR_EST, TOL_LOCAL, DEQPL_SUB
DOUBLE PRECISION :: LAMBDA_DAM_SUB, ETA_CURRENT_SUB
DOUBLE PRECISION :: STRESS_INCREMENT(NTENS), FLOW(NTENS)
DOUBLE PRECISION :: DEV_STRESS(NTENS), DEV_STRESS_NR(NTENS)
DOUBLE PRECISION :: SMEAN, p, p_u, p_u_current, SYIELD_CURRENT
DOUBLE PRECISION :: EE, ENU, FC, FT, A0Y, A1Y, A2Y, A0M, A1M, A2M
DOUBLE PRECISION :: A0R, A1R, A2R, B1, B2, RF, EG, ELAM ! 修复：EG和ELAM为主程序变量
DOUBLE PRECISION :: J2, SMISES, PHI, PHI_NEW, J2_NR, SMISES_NR
DOUBLE PRECISION :: LAMBDA_INC, PC, EPSILON_V_YIELD, DENOMINATOR
DOUBLE PRECISION :: DENOM_Y, DELTA_Y, DENOM_M, DELTA_M
DOUBLE PRECISION :: DENOM_R, DELTA_R, DELTA_SIGMA, dDEQPL
DOUBLE PRECISION :: LAMBDA_M, LAMBDA_R, B3, R_THETA_NR, H
DOUBLE PRECISION :: J3, THETA, PSI, R_PRIME_VAL, ALPHA, PHI_TRIAL
LOGICAL :: IN_TENSILE_SOFTENING, OUTPUT_NEEDED, PLASTIC_FLAG
DOUBLE PRECISION, ALLOCATABLE :: LAMBDA_TABLE(:), ETA_TABLE(:)
DOUBLE PRECISION :: dETA_dLAMBDA, J2_VOL, EPSILON_V, dDEQPL_inc
DOUBLE PRECISION :: Y_eta, H_sigma, B_val, DELTA_M_Y, DELTA_R_M
DOUBLE PRECISION :: DDSDDE_ELASTIC(NTENS, NTENS) ! 弹性刚度矩阵
C CONSTANTS
DOUBLE PRECISION, PARAMETER :: J2_TOL = 1.0D-12,
1 ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0,
2 SIX=6.D0, TOLER=1.0D-6, HALF=0.5D0, THIRD=1.D0/3.D0,
3 TOL_NR_DEF = 1.0D-8
C MATERIAL PARAMETERS
EE = PROPS(1) ! 弹性模量
ENU = PROPS(2) ! 泊松比
FC = PROPS(3) ! 抗压强度
FT = PROPS(4) ! 抗拉强度
A0Y = PROPS(5) ! 初始屈服面参数
A1Y = PROPS(6)
A2Y = PROPS(7)
A0M = PROPS(8) ! 最大屈服面参数
A1M = PROPS(9)
A2M = PROPS(10)
A0R = 0.0D0 ! 残余阶段参数 a0
A1R = PROPS(11) ! 残余阶段参数 a1
A2R = PROPS(12) ! 残余阶段参数 a2
B1 = PROPS(13) ! 压缩损伤参数
B2 = PROPS(14) ! 拉伸损伤参数
RF = PROPS(15) ! 应变率因子
N_ETA_PAIRS = INT(PROPS(16)) ! η-λ表格对数
! 检查KINC异常（修复：避免KINC=0）
IF (KINC <= 0) THEN
WRITE(6,*) "警告：KINC=", KINC, "异常，强制设为1"
KINC = 1
ENDIF
! 输出控制
OUTPUT_NEEDED = .FALSE.
IF (KINC == 1 .OR. MOD(KINC, 100) == 0) THEN
OUTPUT_NEEDED = .TRUE.
ENDIF
! 检查PROPS数组长度
IF (17 + 2*N_ETA_PAIRS + 2 > NPROPS) THEN
WRITE(6,*) "错误：PROPS数组长度不足，至少需要",
1 17 + 2*N_ETA_PAIRS + 2, "个参数"
CALL XIT
ENDIF
! 检查状态变量数组长度
IF (NSTATV < 2*NTENS + 5) THEN
WRITE(6,*) "错误：NSTATV太小，至少需要", 2*NTENS+5
CALL XIT
ENDIF
! 动态分配η-λ表格内存
ALLOCATE(LAMBDA_TABLE(N_ETA_PAIRS), ETA_TABLE(N_ETA_PAIRS))
DO I_ETA = 1, N_ETA_PAIRS
LAMBDA_TABLE(I_ETA) = PROPS(17 + 2*(I_ETA-1)) ! λ1-λ11
ETA_TABLE(I_ETA) = PROPS(18 + 2*(I_ETA-1)) ! η1-η11
ENDDO
LAMBDA_M = PROPS(39) ! λmax
LAMBDA_R = PROPS(40) ! λres
B3 = PROPS(41) ! b3
! 计算弹性刚度矩阵并获取EG和ELAM（修复：传递EG和ELAM到主程序）
CALL INITIAL_STIFFNESS(DDSDDE_ELASTIC, EE, ENU, NTENS, NDI, EG, ELAM)
! 初始化局部变量
STRESS_SUB = 0.0D0
EPLAS_SUB = 0.0D0
LAMBDA_DAM_SUB = 0.0D0
ETA_CURRENT_SUB = 0.0D0
PC = -FT
EPSILON_V_YIELD = 0.0D0
! 检查是否为初始步
IF (TIME(1) <= TOLER) THEN
! 初始步：完全初始化
STATEV(1:NTENS) = 0.0D0 ! 弹性应变
STATEV(NTENS+1:2*NTENS) = 0.0D0 ! 塑性应变
STATEV(2*NTENS+1) = 0.0D0 ! Lambda（损伤累积量）
! 插值获取初始ETA值
CALL INTERPOLATE_ETA(0.0D0, LAMBDA_TABLE,
& ETA_TABLE, N_ETA_PAIRS, ETA_CURRENT_SUB, dETA_dLAMBDA)
STATEV(2*NTENS+2) = ETA_CURRENT_SUB
STATEV(2*NTENS+3) = PC ! 压力截止值
STATEV(2*NTENS+4) = 0.0D0 ! EPSILON_V_YIELD
STATEV(2*NTENS+5) = 0.0D0 ! 历史压力
! 初始应力计算
DO I = 1, NTENS
STRESS(I) = 0.0D0
DO J = 1, NTENS
STRESS(I) = STRESS(I) + DDSDDE_ELASTIC(I, J) * STRAN(J)
ENDDO
ENDDO
STRESS_SUB = STRESS ! 设置子步应力
IF (OUTPUT_NEEDED) THEN
WRITE(6,*) "初始应力 STRESS =", STRESS
WRITE(6,*) "当前步：弹性阶段，无损伤"
WRITE(6,*) "========== 模型初始化完成 =========="
WRITE(6,*) "材料参数: E=", EE, " NU=", ENU, " FC=", FC, " FT=", FT
WRITE(6,*) "初始屈服面参数: A0Y=", A0Y, " A1Y=", A1Y, " A2Y=", A2Y
WRITE(6,*) "剪切模量EG=", EG, " 体积模量ELAM=", ELAM ! 输出验证EG
WRITE(6,*) "--------------------------------------"
ENDIF
ELSE
! 非初始步：从状态变量加载历史数据
DO I = 1, NTENS
EPLAS_SUB(I) = STATEV(NTENS+I) ! 塑性应变
ENDDO
LAMBDA_DAM_SUB = STATEV(2*NTENS+1) ! 损伤累积量
ETA_CURRENT_SUB = STATEV(2*NTENS+2) ! 当前η值
PC = STATEV(2*NTENS+3) ! 压力截止值
EPSILON_V_YIELD = STATEV(2*NTENS+4) ! 体积屈服应变
! STRESS_SUB 从传入的STRESS初始化（上一增量步结束时的应力）
STRESS_SUB = STRESS
IF (OUTPUT_NEEDED) THEN
WRITE(6,*) "从状态变量加载历史数据:"
WRITE(6,*) "塑性应变 =", EPLAS_SUB
WRITE(6,*) "损伤累积量 =", LAMBDA_DAM_SUB
WRITE(6,*) "当前η值 =", ETA_CURRENT_SUB
WRITE(6,*) "当前增量步应变增量 DSTRAN =", DSTRAN
ENDIF
ENDIF
C ==================================================================
C 自适应显式积分主循环
C ==================================================================
! 初始化自适应积分参数
DT_FACTOR = 0.1D0 ! 初始步长因子
TOTAL_STRAIN = 0.0D0
MAX_PLASTIC_INCR = 0.01D0 ! 最大允许塑性应变增量
TOL_LOCAL = 1.0D-6 ! 局部误差容差
N_SUBSTEPS = 0
STRESS_SUB = STRESS
ITER = 0 ! 迭代计数器初始化
PLASTIC_FLAG = .FALSE. ! 标记是否发生塑性变形
! 计算体积应变
EPSILON_V = STRAN(1) + STRAN(2) + STRAN(3) +
& DSTRAN(1) + DSTRAN(2) + DSTRAN(3)
! 主循环：处理整个应变增量
DO WHILE (TOTAL_STRAIN < 1.0D0 .AND. N_SUBSTEPS < 1000)
N_SUBSTEPS = N_SUBSTEPS + 1
ITER = ITER + 1 ! 迭代计数
! 计算当前子步的应变增量
DSTRAN_SUB = DSTRAN * DT_FACTOR
! 弹性预测 - 使用弹性刚度矩阵
DO I = 1, NTENS
STRESS_INCREMENT(I) = 0.0D0
DO J = 1, NTENS
STRESS_INCREMENT(I) = STRESS_INCREMENT(I)
1 + DDSDDE_ELASTIC(I, J) * DSTRAN_SUB(J)
ENDDO
ENDDO
STRESS_TRIAL_SUB = STRESS_SUB + STRESS_INCREMENT
IF (OUTPUT_NEEDED) THEN
WRITE(6,*) "子步", N_SUBSTEPS,
1 "弹性预测应力 STRESS_TRIAL_SUB =", STRESS_TRIAL_SUB
WRITE(6,*) "应变增量 DSTRAN_SUB =", DSTRAN_SUB
ENDIF
! 静水压力计算
SMEAN = (STRESS_TRIAL_SUB(1)+STRESS_TRIAL_SUB(2)+
& STRESS_TRIAL_SUB(3))/THREE
p = -SMEAN ! p > 0为压缩，p < 0为拉伸
! 插值获取当前损伤参数
CALL INTERPOLATE_ETA(LAMBDA_DAM_SUB, LAMBDA_TABLE,
& ETA_TABLE, N_ETA_PAIRS, ETA_CURRENT_SUB, dETA_dLAMBDA)
! 计算屈服应力
IF (RF > 1.0D0 + TOLER) THEN
p_u = p / RF
ELSE
p_u = p
ENDIF
p_u_current = p_u ! 保存用于输出
IN_TENSILE_SOFTENING = (p < 0.0D0) .AND.
& (LAMBDA_DAM_SUB > LAMBDA_M)
IF (IN_TENSILE_SOFTENING) THEN
SYIELD_CURRENT = 3.0D0 * (p + ETA_CURRENT_SUB * RF * FT)
! 修复：拉伸软化区R_THETA_NR赋值（单轴拉伸简化）
R_THETA_NR = 0.5D0
ELSE
DENOM_Y = A1Y + A2Y * p_u
IF (ABS(DENOM_Y) < TOLER) DENOM_Y = SIGN(TOLER, DENOM_Y)
DELTA_Y = A0Y + p_u / DENOM_Y
DENOM_M = A1M + A2M * p_u
DELTA_M = A0M + p_u / DENOM_M
DENOM_R = A1R + A2R * p_u
DELTA_R = A0R + p_u / DENOM_R
IF (LAMBDA_DAM_SUB <= LAMBDA_M) THEN
DELTA_SIGMA = DELTA_Y + ETA_CURRENT_SUB * (DELTA_M - DELTA_Y)
DELTA_M_Y = DELTA_M - DELTA_Y ! 用于计算H
ELSEIF (LAMBDA_DAM_SUB <= LAMBDA_R) THEN
DELTA_SIGMA = DELTA_M + ETA_CURRENT_SUB * (DELTA_R - DELTA_M)
DELTA_R_M = DELTA_R - DELTA_M ! 用于计算H
ELSE
DELTA_SIGMA = DELTA_R
ENDIF
! 计算三参数Willam-Warnke准则参数
CALL CALCULATE_INVARIANTS(STRESS_TRIAL_SUB, NTENS,
& NDI, J2, J3, DEV_STRESS)
THETA = LODE_ANGLE(DEV_STRESS, J2, J3)
PSI = PSI_FUNCTION(p, FT, FC, A0M, A1M, A2M)
R_PRIME_VAL = R_PRIME(PSI, THETA)
SYIELD_CURRENT = RF * DELTA_SIGMA * R_PRIME_VAL
R_THETA_NR = R_PRIME_VAL ! 非拉伸区赋值
ENDIF
! 偏应力计算
DEV_STRESS(1:3) = STRESS_TRIAL_SUB(1:3) - SMEAN
IF (NTENS > 3) THEN
DEV_STRESS(4:NTENS) = STRESS_TRIAL_SUB(4:NTENS)
ENDIF
DEV_STRESS_NR = DEV_STRESS
! 应力不变量计算
J2 = 0.5D0 * (DEV_STRESS(1)**2 + DEV_STRESS(2)**2
& + DEV_STRESS(3)**2)
IF (NTENS > 3) THEN
DO I = 4, NTENS
J2 = J2 + DEV_STRESS(I)**2
ENDDO
ENDIF
J2_NR = J2
! 等效应力
IF (J2 > J2_TOL) THEN
SMISES = SQRT(3.0D0 * J2)
ELSE
SMISES = 0.0D0
ENDIF
SMISES_NR = SMISES
! 屈服判断
PHI = SMISES - SYIELD_CURRENT
PHI_NEW = PHI
! 塑性阶段
IF (PHI.GT.TOLER) THEN
PLASTIC_FLAG = .TRUE.
! 计算硬化模量H
IF (LAMBDA_DAM_SUB <= LAMBDA_M) THEN
H = RF * DELTA_M_Y * dETA_dLAMBDA * R_PRIME_VAL
ELSEIF (LAMBDA_DAM_SUB <= LAMBDA_R) THEN
H = RF * DELTA_R_M * dETA_dLAMBDA * R_PRIME_VAL
ELSE
H = ZERO
ENDIF
! 输出迭代细节
IF (OUTPUT_NEEDED .OR. MOD(ITER, 10) == 0) THEN
WRITE(6,*) "----- 迭代", ITER, "屈服面与硬化模量 -----"
WRITE(6,*) "p_u_current =", p_u_current,
1 " SYIELD_CURRENT =", SYIELD_CURRENT
WRITE(6,*) "H =", H, " EG =", EG, " R_THETA_NR =", R_THETA_NR
WRITE(6,*) "----- 迭代", ITER, "应力修正细节 -----"
WRITE(6,*) "J2_NR =", J2_NR, " SMISES_NR =", SMISES_NR
WRITE(6,*) "PHI_NEW =", PHI_NEW, "目标容差=", TOL_NR_DEF
ENDIF
! 计算流动方向
IF (J2 > J2_TOL) THEN
DO I = 1, NTENS
FLOW(I) = DEV_STRESS(I) / SQRT(2.0D0*J2)
ENDDO
ELSE
FLOW = 0.0D0
ENDIF
! 计算塑性乘子增量（修复：使用正确的EG和R_THETA_NR）
DENOMINATOR = 2.0D0 * EG * R_THETA_NR**2 - H
dDEQPL = PHI_NEW / DENOMINATOR
dDEQPL_inc = DT_FACTOR * dDEQPL
! 输出塑性乘子信息
IF (OUTPUT_NEEDED .OR. MOD(ITER, 10) == 0) THEN
WRITE(6,*) "塑性乘子增量: PHI_NEW=", PHI_NEW,
1 " DENOMINATOR=", DENOMINATOR, " dDEQPL=", dDEQPL
ENDIF
! 更新应力和塑性应变
DEQPL_SUB = ABS(dDEQPL)
STRESS_SUB = STRESS_TRIAL_SUB - 2.0D0 * EG * DEQPL_SUB * FLOW
EPLAS_SUB = EPLAS_SUB + DEQPL_SUB * FLOW
! 计算损伤增量
J2_VOL = 0.5D0 * (DEV_STRESS(1)**2 +
& DEV_STRESS(2)**2 + DEV_STRESS(3)**2)
CALL DamageEvolution(B1, B2, B3, FC, FT, SMEAN,
& DEQPL_SUB, NTENS, LAMBDA_INC, J2_VOL, EPSILON_V,
& EPSILON_V_YIELD, LAMBDA_DAM_SUB, LAMBDA_M, RF)
! 更新损伤参数
LAMBDA_DAM_SUB = LAMBDA_DAM_SUB + LAMBDA_INC
! 修复：更新体积屈服应变（当前体积应变作为新阈值）
EPSILON_V_YIELD = EPSILON_V
ELSE
! 弹性阶段
IF (OUTPUT_NEEDED .AND. N_SUBSTEPS == 1) THEN
WRITE(6,*) "当前步：弹性阶段，使用弹性刚度矩阵"
ENDIF
ENDIF
! 误差估计和步长调整
ERROR_EST = DEQPL_SUB
IF (ERROR_EST < TOL_LOCAL * 0.1D0 .AND. ERROR_EST > ZERO) THEN
DT_FACTOR = MIN(DT_FACTOR * 1.5D0, 1.0D0)
ELSE IF (ERROR_EST > TOL_LOCAL) THEN
DT_FACTOR = MAX(DT_FACTOR * 0.8D0, 1.0D-5)
ENDIF
TOTAL_STRAIN = TOTAL_STRAIN + DT_FACTOR
IF (TOTAL_STRAIN + DT_FACTOR > 1.0D0) THEN
DT_FACTOR = 1.0D0 - TOTAL_STRAIN
ENDIF
ENDDO
C ==================================================================
C 结束自适应积分主循环
C ==================================================================
! 更新最终状态变量
STRESS = STRESS_SUB
STATEV(1:NTENS) = STRAN + DSTRAN - EPLAS_SUB ! 弹性应变
STATEV(NTENS+1:2*NTENS) = EPLAS_SUB ! 塑性应变
STATEV(2*NTENS+1) = LAMBDA_DAM_SUB ! 损伤累积量
CALL INTERPOLATE_ETA(LAMBDA_DAM_SUB, LAMBDA_TABLE,
& ETA_TABLE, N_ETA_PAIRS, ETA_CURRENT_SUB, dETA_dLAMBDA)
STATEV(2*NTENS+2) = ETA_CURRENT_SUB
SMEAN = (STRESS(1)+STRESS(2)+STRESS(3))/THREE
p = -SMEAN
IF (LAMBDA_DAM_SUB > LAMBDA_M .AND. p < 0.0D0) THEN
PC = -ETA_CURRENT_SUB * FT
STATEV(2*NTENS+3) = PC
ENDIF
STATEV(2*NTENS+4) = EPSILON_V_YIELD ! 保存更新后的体积屈服应变
! 更新刚度矩阵
IF (PLASTIC_FLAG) THEN
CALL CONSISTENT_TANGENT(DDSDDE_ELASTIC, FLOW, EG, H,
& NTENS, NDI, DDSDDE)
IF (OUTPUT_NEEDED) THEN
WRITE(6,*) "使用一致性切线刚度矩阵"
ENDIF
ELSE
DDSDDE = DDSDDE_ELASTIC
IF (OUTPUT_NEEDED) THEN
WRITE(6,*) "使用弹性刚度矩阵"
ENDIF
ENDIF
! 释放内存
DEALLOCATE(LAMBDA_TABLE, ETA_TABLE)
! 增量步结束信息
IF (OUTPUT_NEEDED) THEN
WRITE(6,*) "========== 增量步: ", KINC, " 结束 =========="
WRITE(6,*) "总子步数: ", N_SUBSTEPS, " 最终迭代次数: ", ITER
WRITE(6,*) "最终损伤值: ", LAMBDA_DAM_SUB, " 最终等效应力: ", SMISES_NR
WRITE(6,*) "--------------------------------------"
ENDIF
C-----------------------------------------------------------------------
C 内部子程序定义
C-----------------------------------------------------------------------
CONTAINS
SUBROUTINE INITIAL_STIFFNESS(DDSDDE, EE, ENU, NTENS, NDI, EG, ELAM)
! 修复：将EG和ELAM作为输出参数，传递到主程序
IMPLICIT NONE
INTEGER, INTENT(IN) :: NTENS, NDI
DOUBLE PRECISION, INTENT(IN) :: EE, ENU
DOUBLE PRECISION, INTENT(OUT) :: DDSDDE(NTENS, NTENS), EG, ELAM
INTEGER :: I, J
EG = EE / (2.0D0*(1.0D0 + ENU)) ! 剪切模量计算
ELAM = EE*ENU / ((1.0D0 + ENU)*(1.0D0 - 2.0D0*ENU)) ! 体积模量
DDSDDE = 0.0D0
DO I = 1, NDI
DO J = 1, NDI
DDSDDE(I,J) = ELAM
END DO
DDSDDE(I,I) = ELAM + 2.0D0*EG
END DO
DO I = NDI+1, NTENS
DDSDDE(I,I) = EG
END DO
END SUBROUTINE INITIAL_STIFFNESS
SUBROUTINE DamageEvolution(B1, B2, B3, FC, FT, SMEAN_NR,
& dDEQPL_inc, NTENS, LAMBDA_INC, J2_VOL, EPSILON_V,
& EPSILON_V_YIELD, LAMBDA_CURRENT, LAMBDA_M, RF)
IMPLICIT NONE
INTEGER, INTENT(IN) :: NTENS
DOUBLE PRECISION, INTENT(IN) :: B1, B2, FC, FT, LAMBDA_CURRENT,
1 SMEAN_NR, dDEQPL_inc, J2_VOL, EPSILON_V, B3,
2 EPSILON_V_YIELD, LAMBDA_M, RF
DOUBLE PRECISION :: P_NR, RATIO, FD, denominator, temp
DOUBLE PRECISION, INTENT(OUT) :: LAMBDA_INC
DOUBLE PRECISION :: LAMBDA_INC_VOL = 0.0D0
P_NR = -SMEAN_NR
LAMBDA_INC = 0.0D0
IF (P_NR < 0.0D0) THEN ! 拉伸区（p < 0）
temp = 1.0D0 - ABS(P_NR)/(RF*FT)
temp = MAX(temp, 1.0D-3)
denominator = RF * temp**B2
ELSE ! 压缩区（p ≥ 0）
temp = 1.0D0 + P_NR/(RF*FC)
denominator = RF * temp**B1
ENDIF
IF (denominator < 1.0D-20) THEN
denominator = 1.0D-10
ENDIF
LAMBDA_INC = dDEQPL_inc / denominator
! 体积损伤项
IF (J2_VOL < 1.0D-12 .AND. P_NR < 0.0D0) THEN
RATIO = SQRT(3.0D0*J2_VOL) / (ABS(P_NR) + 1.0D-12)
IF (RATIO < 0.1D0) THEN
FD = 1.0D0 - RATIO/0.1D0
LAMBDA_INC_VOL = B3 * FD * (EPSILON_V - EPSILON_V_YIELD)
LAMBDA_INC = LAMBDA_INC + LAMBDA_INC_VOL
ENDIF
ENDIF
END SUBROUTINE DamageEvolution
SUBROUTINE INTERPOLATE_ETA(LAMBDA_CURRENT,
1 LAMBDA_TABLE, ETA_TABLE, N, ETA, dETA_dLAMBDA)
IMPLICIT NONE
INTEGER, INTENT(IN) :: N
DOUBLE PRECISION, INTENT(IN) :: LAMBDA_CURRENT,
1 LAMBDA_TABLE(N), ETA_TABLE(N)
DOUBLE PRECISION, INTENT(OUT) :: ETA, dETA_dLAMBDA
INTEGER :: I
! 修复：初始化dETA_dLAMBDA为0（处理N=1的情况）
dETA_dLAMBDA = 0.0D0
IF (LAMBDA_CURRENT <= LAMBDA_TABLE(1)) THEN
ETA = ETA_TABLE(1)
IF (N >= 2) THEN
dETA_dLAMBDA = (ETA_TABLE(2) - ETA_TABLE(1)) /
& (LAMBDA_TABLE(2) - LAMBDA_TABLE(1))
ENDIF
RETURN
ENDIF
IF (LAMBDA_CURRENT >= LAMBDA_TABLE(N)) THEN
ETA = ETA_TABLE(N)
IF (N >= 2) THEN
dETA_dLAMBDA = (ETA_TABLE(N) - ETA_TABLE(N-1)) /
& (LAMBDA_TABLE(N) - LAMBDA_TABLE(N-1))
ENDIF
RETURN
ENDIF
DO I = 1, N-1
IF (LAMBDA_CURRENT >= LAMBDA_TABLE(I) .AND.
& LAMBDA_CURRENT <= LAMBDA_TABLE(I+1)) THEN
ETA = ETA_TABLE(I) + (LAMBDA_CURRENT - LAMBDA_TABLE(I)) *
& (ETA_TABLE(I+1) - ETA_TABLE(I)) /
& (LAMBDA_TABLE(I+1) - LAMBDA_TABLE(I))
dETA_dLAMBDA = (ETA_TABLE(I+1) - ETA_TABLE(I)) /
& (LAMBDA_TABLE(I+1) - LAMBDA_TABLE(I))
RETURN
ENDIF
ENDDO
ETA = ETA_TABLE(1)
dETA_dLAMBDA = 0.0D0
END SUBROUTINE INTERPOLATE_ETA
SUBROUTINE CALCULATE_INVARIANTS(STRESS, NTENS, NDI, J2, J3, DEV)
IMPLICIT NONE
INTEGER, INTENT(IN) :: NTENS, NDI
DOUBLE PRECISION, INTENT(IN) :: STRESS(NTENS)
DOUBLE PRECISION, INTENT(OUT) :: J2, J3, DEV(NTENS)
DOUBLE PRECISION :: SMEAN
INTEGER :: I
SMEAN = (STRESS(1) + STRESS(2) + STRESS(3)) / 3.0D0
DEV(1:3) = STRESS(1:3) - SMEAN
IF (NTENS > 3) THEN
DEV(4:NTENS) = STRESS(4:NTENS)
ENDIF
J2 = 0.5D0 * (DEV(1)**2 + DEV(2)**2 + DEV(3)**2)
IF (NTENS > 3) THEN
DO I = 4, NTENS
J2 = J2 + DEV(I)**2
ENDDO
ENDIF
J3 = DEV(1)*DEV(2)*DEV(3)
IF (NTENS >= 6) THEN
J3 = J3 - DEV(1)*DEV(6)**2 - DEV(2)*DEV(5)**2 - DEV(3)*DEV(4)**2
1 + 2.0D0*DEV(4)*DEV(5)*DEV(6)
ENDIF
END SUBROUTINE CALCULATE_INVARIANTS
DOUBLE PRECISION FUNCTION R_PRIME(PSI, THETA)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: PSI, THETA
DOUBLE PRECISION :: COST, NUM, DEN
COST = COS(THETA)
NUM = 2.0D0*(1.0D0 - PSI**2)*COST + (2.0D0*PSI - 1.0D0) *
1 SQRT(4.0D0*(1.0D0 - PSI**2)*COST**2 + 5.0D0*PSI**2 - 4.0D0*PSI)
DEN = 4.0D0*(1.0D0 - PSI**2)*COST**2 + (1.0D0 - 2.0D0*PSI)**2
R_PRIME = NUM / DEN
END FUNCTION R_PRIME
DOUBLE PRECISION FUNCTION PSI_FUNCTION(p, FT, FC, A0M, A1M, A2M)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: p, FT, FC, A0M, A1M, A2M
DOUBLE PRECISION :: ALPHA = 1.15D0, PSI_TEMP
IF (p <= 0.0D0) THEN
PSI_TEMP = 0.5D0
ELSE IF (p <= FC/3.0D0) THEN
PSI_TEMP = 0.5D0 + 1.5D0 * FT / FC
ELSE IF (p <= 2.0D0*ALPHA*FC/3.0D0) THEN
PSI_TEMP = ALPHA * FC / (A0M + (2.0D0*ALPHA*FC/3.0D0)
1 / (A1M + A2M*2.0D0*ALPHA*FC/3.0D0))
ELSE IF (p <= 3.0D0*FC) THEN
PSI_TEMP = 0.753D0
ELSE
PSI_TEMP = 1.0D0
END IF
! 限制PSI范围，避免异常
PSI_TEMP = MAX(MIN(PSI_TEMP, 1.2D0), 0.1D0)
PSI_FUNCTION = PSI_TEMP
END FUNCTION PSI_FUNCTION
DOUBLE PRECISION FUNCTION LODE_ANGLE(DEV_STRESS, J2, J3)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: DEV_STRESS(6), J2, J3
DOUBLE PRECISION :: COS3THETA
IF (J2 < 1.0D-12) THEN
LODE_ANGLE = 0.0D0
RETURN
END IF
COS3THETA = (3.0D0 * SQRT(3.0D0) * J3) / (2.0D0 * J2**(1.5D0))
COS3THETA = MAX(MIN(COS3THETA, 1.0D0), -1.0D0)
LODE_ANGLE = ACOS(COS3THETA) / 3.0D0
END FUNCTION LODE_ANGLE
SUBROUTINE CONSISTENT_TANGENT(DDSDDE_ELASTIC, FLOW, EG, H,
& NTENS, NDI, DDSDDE_EP)
IMPLICIT NONE
INTEGER, INTENT(IN) :: NTENS, NDI
DOUBLE PRECISION, INTENT(IN) :: DDSDDE_ELASTIC(NTENS, NTENS)
DOUBLE PRECISION, INTENT(IN) :: FLOW(NTENS), EG, H
DOUBLE PRECISION, INTENT(OUT) :: DDSDDE_EP(NTENS, NTENS)
INTEGER :: I, J
DOUBLE PRECISION :: denom, scale_factor
DOUBLE PRECISION, PARAMETER :: ZERO=0.D0, ONE=1.D0, TWO=2.D0,
& THREE=3.D0, TOL=1.D-12
DDSDDE_EP = DDSDDE_ELASTIC
denom = 2.0D0 * EG - H
IF (ABS(denom) < TOL) THEN
denom = SIGN(TOL, denom)
END IF
scale_factor = (4.0D0 * EG * EG) / denom
DO I = 1, NTENS
DO J = 1, NTENS
DDSDDE_EP(I, J) = DDSDDE_EP(I, J) -
& scale_factor * FLOW(I) * FLOW(J)
END DO
END DO
END SUBROUTINE CONSISTENT_TANGENT
END SUBROUTINE UMAT
扫一扫
1360
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
