3. 弹性本构UMAT实现

各向同性线弹性是UMAT开发中最基础、也最适合用来打通“公式—代码—验证”闭环的本构模型。它的物理含义清晰：应力与弹性应变之间通过四阶弹性张量线性映射，材料只有两个独立常数，通常取弹性模量 E 和泊松比 ν，或取体积模量 K 与剪切模量 G。本节从应力更新公式出发，推导Jacobian矩阵，给出可直接编译的Fortran实现，并用单轴拉伸算例验证。

一、应力更新公式

在增量有限元中，UMAT在每个积分点接收上一步收敛的应力 STRESS、总应变增量 DSTRAN，以及状态变量 STATEV。线弹性无历史依赖，因此应力更新为：

σ(n+1) = σ(n) + C : Δε

其中 C 为四阶弹性张量。对三维问题，用Lamé常数表示最简洁：

λ = Eν / [(1+ν)(1-2ν)]
μ = G = E / [2(1+ν)]
C_ijkl = λ δ_ij δ_kl + μ (δ_ik δ_jl + δ_il δ_jk)

展开为Voigt记法的6×6矩阵（工程剪应变约定）：

D11 = D22 = D33 = λ + 2μ
D12 = D13 = D23 = λ
D44 = D55 = D66 = μ

其余分量为零。注意：ABAQUS的UMAT中剪应变分量为工程剪应变 γ = 2ε，因此剪切刚度项直接为 μ，不需要额外乘2。这是初学者最容易出错的地方之一。

二、Jacobian矩阵推导

Jacobian矩阵 DDSDDE(i,j) 定义为应力增量对应变增量的偏导：

DDSDDE(i,j) = ∂Δσ_i / ∂Δε_j

对线弹性本构，应力增量与应变增量严格线性，因此Jacobian就是弹性张量本身，且在整个增量步内为常数。这意味着：

  • Jacobian不随当前应力或应变变化，无需迭代更新；
  • 它是对称正定矩阵，保证Newton迭代的收敛性；
  • 在UMAT中只需在增量步开始时赋值一次。

若采用 K-G 形式，矩阵元素为：

D11 = K + 4G/3
D12 = K - 2G/3
D44 = G

两种写法等价，选择哪种取决于输入参数。工程中常直接给 E、ν，因此代码里先算 λ、μ 更自然。

各向同性线弹性UMAT：输入—更新—输出与Jacobian结构 UMAT输入 PROPS(1)=E PROPS(2)=ν STRESS(6)：σ(n) DSTRAN(6)：Δε STATEV(*)：本构无依赖 中间计算 λ = Eν/[(1+ν)(1-2ν)] μ = E/[2(1+ν)] D11 = λ+2μ D12 = λ D44 = μ UMAT输出 STRESS(6)：σ(n+1) DDSDDE(6,6) STATEV：可留空 SSE, SPD, SCD DDSDDE 矩阵结构（Voigt记法，对称） ε11 ε22 ε33 γ12 γ13 γ23 σ11 σ22 σ33 τ12 τ13 λ+2μ 与 λ 填充 μ 填充 0 0 线弹性：DDSDDE为常数，与σ、ε无关；剪应变采用工程剪应变约定

三、Fortran代码实现

下面给出一个完整的UMAT子程序，适用于三维实体单元。代码中所有变量按ABAQUS约定声明，DDSDDE 先清零再赋值。

      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)
      INCLUDE 'ABA_PARAM.INC'
      CHARACTER*80 CMNAME
      DIMENSION 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)

      REAL*8 E, ENU, ELAM, EMU, D11, D12

      E   = PROPS(1)
      ENU = PROPS(2)

      ELAM = E*ENU/((1.0D0+ENU)*(1.0D0-2.0D0*ENU))
      EMU  = E/(2.0D0*(1.0D0+ENU))
      D11  = ELAM + 2.0D0*EMU
      D12  = ELAM

      DO I=1,NTENS
        DO J=1,NTENS
          DDSDDE(I,J) = 0.0D0
        END DO
      END DO

      DO I=1,3
        DO J=1,3
          IF (I.EQ.J) THEN
            DDSDDE(I,J) = D11
          ELSE
            DDSDDE(I,J) = D12
          END IF
        END DO
      END DO

      DO I=4,NTENS
        DDSDDE(I,I) = EMU
      END DO

      DO I=1,NTENS
        DO J=1,NTENS
          STRESS(I) = STRESS(I) + DDSDDE(I,J)*DSTRAN(J)
        END DO
      END DO

      SSE = 0.0D0
      DO I=1,NTENS
        SSE = SSE + 0.5D0*STRESS(I)*STRAN(I)
      END DO
      SPD = 0.0D0
      SCD = 0.0D0

      RETURN
      END

代码要点：NTENS 对三维单元为6，对平面应变单元为4，对平面应力单元为3。若用于平面应力，需将 DDSDDE 改为平面应力矩阵，或采用ABAQUS的平面应力单元自动处理。本代码按三维/平面应变通用写法，NDI=3、NSHR=3 时正确。

四、单轴拉伸验证

验证模型：一个100 mm×10 mm×10 mm的立方体单元，左端固定，右端施加 u_x = 0.1 mm 的位移。材料参数取 E = 210000 MPa、ν = 0.3。理论解为：

ε_x = 0.1/100 = 1.0e-3
σ_x = E ε_x = 210 MPa
ε_y = ε_z = -ν ε_x = -3.0e-4

在ABAQUS中提交后，检查积分点的 STRESS(1) 是否等于210 MPa，STRESS(2)、STRESS(3) 是否接近0（自由边界），STRAN(2)、STRAN(3) 是否为-3.0e-4。若偏差超过1%，优先检查单位一致性、泊松比符号、以及是否误用了工程剪应变。

验证清单：
  • 轴向应力 σ_x = E·ε_x，误差应小于0.1%；
  • 横向应变 ε_y = ε_z = -ν·ε_x；
  • DDSDDE 对称且正定，行列式大于0；
  • SSE 与 0.5·σ:ε 一致，量级合理；
  • 改变网格密度，结果不变（线弹性无网格依赖）。

通过单轴拉伸后，可进一步做纯剪切验证：施加 γ_xy = 1.0e-3，检查 τ_xy = G·γ_xy，其中 G = E/[2(1+ν)] = 80769.2 MPa。这一步能有效暴露剪应变约定错误。完成这两个算例，各向同性线弹性UMAT即可视为正确，可作为后续塑性、损伤等本构的开发基线。

← 上一章 📖 返回目录 下一章 →