01 概述 / Overview Overview

1.1 本构模型简介 / Constitutive Model Introduction 1.1 Constitutive Model Introduction

Huang（Yonggang Huang）教授开发的UMAT单晶塑性子程序是晶体塑性有限元方法（CPFEM）领域最广泛使用的开源实现之一。 该子程序基于Abaqus/Standard用户材料子程序（UMAT）接口，采用隐式时间积分方案， 能够准确模拟单晶及多晶材料在复杂加载条件下的塑性变形行为。 该模型考虑了晶体的各向异性弹性、晶体学滑移系上的剪切变形，以及位错交互作用导致的硬化效应。

The UMAT single-crystal plasticity subroutine developed by Professor Yonggang Huang is one of the most widely used open-source implementations in the field of Crystal Plasticity Finite Element Method (CPFEM). Based on the Abaqus/Standard User Material (UMAT) interface with an implicit time integration scheme, it accurately simulates the plastic deformation behavior of single-crystal and polycrystalline materials under complex loading conditions. The model accounts for anisotropic elasticity, crystallographic slip system shear deformation, and hardening effects due to dislocation interactions.

本程序的核心特征包括：基于中间构型的乘法分解运动学框架、Schmid定律驱动的滑移系激活判据、 率相关幂律流动法则、可选择的硬化模型（双曲正割硬化或Bassani-Wu硬化）， 以及完整的隐式应力更新算法和一致切线刚度矩阵（Jacobian）计算。

Key features of this subroutine include: multiplicative decomposition kinematics based on the intermediate configuration, Schmid’s law driven slip system activation criterion, rate-dependent power-law flow rule, selectable hardening models (hyperbolic secant hardening or Bassani-Wu hardening), and a complete implicit stress update algorithm with consistent tangent stiffness matrix (Jacobian) computation.

1.2 理论基础 / Theoretical Foundation 1.2 Theoretical Foundation

黄永刚UMAT的理论基础主要来源于Peirce、Asaro和Needleman（1983）以及Peirce、Shih和Needleman（1984）的经典工作。 这些文献建立了晶体塑性理论在现代有限元框架下的数学表述，包括：

The theoretical foundation of Huang’s UMAT is primarily derived from the seminal works of Peirce, Asaro and Needleman (1983) and Peirce, Shih and Needleman (1984). These papers established the mathematical formulation of crystal plasticity theory within the modern finite element framework, including:

• 总变形梯度F的乘法分解为弹性部分Fe和塑性部分Fp Multiplicative decomposition of the total deformation gradient F into elastic Fe and plastic Fp parts
• 基于Schmid分解的分解应力（resolved shear stress）概念 The concept of resolved shear stress based on Schmid decomposition
• 率相关流动法则，避免率无关模型的非唯一性和数值不稳定性 Rate-dependent flow rule to avoid non-uniqueness and numerical instability of rate-independent models
• 自硬化（self-hardening）和潜硬化（latent hardening）的统一描述 Unified description of self-hardening and latent hardening

后续Asaro和Needleman（1985）以及Harren、Deve和Asaro（1988）的工作进一步完善了晶体塑性理论在有限变形条件下的数值实现。 黄永刚在其子程序中整合了这些理论成果，并针对Abaqus/Standard的隐式求解器进行了算法优化。

Subsequent works by Asaro and Needleman (1985) and Harren, Deve and Asaro (1988) further refined the numerical implementation of crystal plasticity theory under finite deformation conditions. Huang integrated these theoretical achievements into his subroutine and optimized the algorithm specifically for the Abaqus/Standard implicit solver.

1.3 适用范围 / Scope 1.3 Scope

本子程序适用于以下研究场景：

This subroutine is applicable to the following research scenarios:

• FCC、BCC和HCP晶体结构的单晶和多晶金属材料的力学行为模拟 Mechanical behavior simulation of single-crystal and polycrystalline metals with FCC, BCC, and HCP crystal structures
• 单调加载、循环加载及应变路径变化条件下的塑性变形分析 Plastic deformation analysis under monotonic, cyclic, and strain-path-change loading conditions
• 织构演化（texture evolution）和塑性各向异性发展的定量预测 Quantitative prediction of texture evolution and plastic anisotropy development
• 局部化变形、剪切带形成及细观尺度损伤起始研究 Studies of localized deformation, shear band formation, and mesoscale damage initiation
• 微柱压缩、纳米压痕等小尺度单晶实验的有限元建模 Finite element modeling of micro-pillar compression, nanoindentation, and other small-scale single-crystal experiments

02 本构理论 / Constitutive Theory Constitutive Theory

2.1 分解运动学 / Kinematics of Decomposition 2.1 Kinematics of Decomposition

晶体塑性理论采用乘法分解（multiplicative decomposition）将总变形梯度F分解为弹性部分Fe和塑性部分Fp：

Crystal plasticity theory employs multiplicative decomposition to separate the total deformation gradient F into elastic Fe and plastic Fp parts:

其中，塑性变形梯度Fp由所有激活滑移系上的剪切应变累积产生，保持晶格结构不变（等容塑性变形）； 弹性变形梯度Fe描述晶格的弹性畸变和刚性转动。 该分解将参考构型（reference configuration）映射到中间构型（intermediate configuration）， 再映射到当前构型（current configuration）。

Here, the plastic deformation gradient Fp is produced by the accumulation of shear strains on all active slip systems, preserving the crystal lattice structure (isochoric plastic deformation); the elastic deformation gradient Fe describes the elastic distortion and rigid rotation of the lattice. This decomposition maps from the reference configuration to the intermediate configuration, and then to the current configuration.

塑性变形梯度的演化方程为：

The evolution equation for the plastic deformation gradient is:

上式中，s^(α)和m^(α)分别是第α个滑移系的滑移方向（slip direction）和滑移面法向（slip plane normal）， 二者定义在中间构型中且保持正交归一。γ̇^(α)为对应的剪切应变率。 滑移方向s和法向m在当前构型中的映像分别为F^es和(F^e)^-Tm。

In the above equation, s^(α) and m^(α) are the slip direction and slip plane normal of the α-th slip system, respectively, defined in the intermediate configuration and maintained orthonormal. γ̇^(α) is the corresponding shear strain rate. The images of slip direction s and normal m in the current configuration are F^es and (F^e)^-Tm, respectively.

弹性变形梯度F^e的极分解为F^e = R^eU^e = V^eR^e， 其中U^e和V^e分别为右、左弹性拉伸张量，R^e为弹性转动张量。 Green弹性应变张量定义为：

The polar decomposition of the elastic deformation gradient is F^e = R^eU^e = V^eR^e, where U^e and V^e are the right and left elastic stretch tensors, and R^e is the elastic rotation tensor. The Green elastic strain tensor is defined as:

弹性本构关系采用广义胡克定律，以第二类Piola-Kirchhoff应力S^e（定义在中间构型）和Green弹性应变E^e表示：

The elastic constitutive relation adopts the generalized Hooke’s law, expressed in terms of the second Piola-Kirchhoff stress S^e (defined on the intermediate configuration) and Green elastic strain E^e:

当前构型中的Kirchhoff应力τ和Cauchy应力σ分别为：

The Kirchhoff stress τ and Cauchy stress σ in the current configuration are:

2.2 Schmid定律 / Schmid’s Law 2.2 Schmid’s Law

Schmid定律是晶体塑性理论的核心，它指出：滑移系α上的分解剪应力（resolved shear stress）τ^(α)等于Cauchy应力σ在滑移方向-滑移面法向投影上的分量。 在程序实现中，更常用的是对称化的Schmid张量P^(α)：

Schmid’s law is the core of crystal plasticity theory, stating that the resolved shear stress τ^(α) on slip system α equals the component of Cauchy stress σ projected onto the slip direction-slip plane normal pair. In program implementation, the symmetrized Schmid tensor P^(α) is more commonly used:

分解剪应力可表示为：

The resolved shear stress can be expressed as:

这里”:”表示双点积（张量缩并），即σ_ijP_ij。Schmid张量P^(α)是对称张量， 它将应力空间中驱动第α个滑移系剪切变形的方向提取出来。 当τ^(α)达到该滑移系的临界分切应力（critical resolved shear stress, CRSS）g^(α)时，滑移开始激活。

Here “:” denotes the double dot product (tensor contraction), i.e., σ_ijP_ij. The Schmid tensor P^(α) is symmetric and extracts the direction in stress space that drives shear deformation on the α-th slip system. Slip activates when τ^(α) reaches the critical resolved shear stress (CRSS) g^(α) of that slip system.

2.3 率相关流动律 / Rate-Dependent Flow Law 2.3 Rate-Dependent Flow Law

黄永刚UMAT采用率相关（viscoplastic）流动法则，避免了率无关模型中滑移系激活-非激活切换带来的数值困难。 每个滑移系上的剪切应变率由以下幂律（power law）描述：

Huang’s UMAT employs a rate-dependent (viscoplastic) flow rule, avoiding numerical difficulties associated with the activation-deactivation switching of slip systems in rate-independent models. The shear strain rate on each slip system is described by the following power law:

式中各参数含义如下：

The parameters in the equation are defined as follows:

• γ̇₀: 参考剪切率（reference shear rate），单位为s^-1，通常取与实验应变速率相同量级 Reference shear rate, in s^-1, typically of the same order as the experimental strain rate
• n: 率敏感指数（rate sensitivity exponent），无量纲。n值越大，材料行为越接近率无关；n→∞时退化为率无关模型 Rate sensitivity exponent, dimensionless. Larger n means behavior closer to rate-independent; n→∞ recovers the rate-independent model
• τ^(α): 第α个滑移系上的分解剪应力 Resolved shear stress on the α-th slip system
• g^(α): 第α个滑移系的当前临界分切应力（即当前强度/硬化状态） Current critical resolved shear stress (i.e., current strength/hardening state) of the α-th slip system
• sign(x): 符号函数，取值为-1、0或1 Sign function, taking values -1, 0, or 1

率相关模型的优势在于：所有滑移系始终处于”部分激活”状态，不需要显式的激活/非激活切换条件； 这大大提高了数值计算的稳定性和收敛性。 在典型金属中，n值范围约为20-200：FCC金属n≈20-100，BCC金属n≈50-200，HCP金属n≈20-80。 n=1对应线性粘性（牛顿流体），n很大时近似于理想塑性。

The advantage of the rate-dependent model is that all slip systems are always in a “partially active” state, without requiring explicit active/inactive switching conditions; this greatly improves numerical stability and convergence. In typical metals, n ranges from about 20-200: FCC metals n≈20-100, BCC metals n≈50-200, HCP metals n≈20-80. n=1 corresponds to linear viscosity (Newtonian fluid), and large n approximates ideal plasticity.

2.4 硬化律 / Hardening Laws 2.4 Hardening Laws

临界分切应力g^(α)随累积塑性剪切应变γ = Σ_α ∫|γ̇^(α)|dt而演化， 演化方程为：

The critical resolved shear stress g^(α) evolves with the accumulated plastic shear strain γ = Σ_α ∫|γ̇^(α)|dt, according to:

其中h_αβ为硬化矩阵，其对角元素h_αα表示自硬化（self-hardening）， 非对角元素h_αβ(α≠β)表示潜硬化（latent hardening）， 描述不同滑移系之间的交互强化效应。 黄永刚UMAT提供两种硬化模型选项：

where h_αβ is the hardening matrix. The diagonal elements h_αα represent self-hardening, and the off-diagonal elements h_αβ (α≠β) represent latent hardening, describing the interaction strengthening effect between different slip systems. Huang’s UMAT provides two hardening model options:

(a) 双曲正割硬化模型 / Hyperbolic Secant Hardening Model

自硬化模量采用双曲正割（sech²）形式：

The self-hardening modulus adopts a hyperbolic secant (sech²) form:

累积剪切应变γ从0增加到∞时，硬化模量h从H₀单调递减到0， 临界分切应力从初始值τ₀渐近趋近于饱和值τ_s。 积分后可得g(γ)的显式表达式：

As the accumulated shear strain γ increases from 0 to ∞, the hardening modulus h decreases monotonically from H₀ to 0, and the critical resolved shear stress asymptotically approaches the saturation value τ_s from the initial value τ₀. Integration yields the explicit expression for g(γ):

参数说明：H₀为初始硬化模量（MPa），τ₀为初始CRSS（MPa），τ_s为饱和CRSS（MPa）。 该模型简洁且物理意义明确，适用于大多数FCC金属。

Parameter definitions: H₀ is the initial hardening modulus (MPa), τ₀ is the initial CRSS (MPa), and τ_s is the saturation CRSS (MPa). This model is concise with clear physical meaning and is suitable for most FCC metals.

(b) Bassani-Wu硬化模型 / Bassani-Wu Hardening Model

Bassani和Wu（1991）提出的硬化模型考虑了滑移系间交互作用的各向异性特征。 自硬化模量F(γ)为：

The hardening model proposed by Bassani and Wu (1991) accounts for the anisotropic characteristics of interaction between slip systems. The self-hardening modulus F(γ) is:

潜硬化系数q_αβ取决于滑移系对(α,β)的交互类型（共面/共向/正交/Lomer-Cottrell锁等）， 硬化矩阵一般写作h_αβ = q_αβ · h_ββ。 完整的潜硬化矩阵是12×12的对称矩阵，包含7种不同的交互类型。

The latent hardening coefficient q_αβ depends on the interaction type of the slip system pair (α,β) (coplanar/codirectional/orthogonal/Lomer-Cottrell lock, etc.). The hardening matrix is generally written as h_αβ = q_αβ · h_ββ. The full latent hardening matrix is a 12×12 symmetric matrix containing 7 different interaction types.

Bassani模型的优势在于能够更准确地描述FCC晶体中不同滑移系之间的交互硬化效应， 特别是非共面滑移系之间的强相互作用（如Lomer-Cottrell锁）。 但该模型需要额外的潜硬化参数，且交互类型分类较为复杂。

The advantage of the Bassani model is its ability to more accurately describe the interaction hardening effects between different slip systems in FCC crystals, especially the strong interactions between non-coplanar slip systems (such as Lomer-Cottrell locks). However, this model requires additional latent hardening parameters, and the interaction type classification is relatively complex.

2.5 应力更新 / Stress Update (Jaumann Rate) 2.5 Stress Update (Jaumann Rate)

在有限变形框架下，本构关系必须以客观（objective）的应力率形式表述，以消除刚性转动带来的虚假应力变化。 黄永刚UMAT采用Jaumann应力率：

In the finite deformation framework, the constitutive relation must be expressed in an objective stress rate form to eliminate spurious stress changes due to rigid body rotation. Huang’s UMAT adopts the Jaumann stress rate:

其中∇σ = σ̇ – W·σ + σ·W为Jaumann应力率， D为变形率（symmetric part of velocity gradient L）， W为旋率（skew-symmetric part of L）， C^ep为弹塑性刚度张量，D^p为塑性变形率：

where ∇σ = σ̇ – W·σ + σ·W is the Jaumann stress rate, D is the rate of deformation (symmetric part of velocity gradient L), W is the spin (skew-symmetric part of L), C^ep is the elastoplastic stiffness tensor, and D^p is the plastic rate of deformation:

隐式时间积分采用Newton-Raphson迭代求解每个时间步末的应力和状态变量。 线性化后得到以Δγ^(α)为未知量的线性方程组，通过LU分解求解。 收敛后更新STATEV数组中的应力、应变和各滑移系变量。

Implicit time integration employs Newton-Raphson iteration to solve for the stress and state variables at the end of each time increment. Linearization yields a linear system of equations with Δγ^(α) as unknowns, solved by LU decomposition. After convergence, the STATEV array is updated with stress, strain, and slip system variables.

03 Abaqus UMAT接口 / UMAT Interface UMAT Interface

3.1 完整Fortran签名 / Full Fortran Signature 3.1 Full Fortran Signature

Abaqus/Standard的用户材料子程序UMAT具有固定的接口签名。黄永刚单晶塑性子程序在该框架下实现，完整声明如下：

The Abaqus/Standard User Material subroutine (UMAT) has a fixed interface signature. Huang’s single-crystal plasticity subroutine is implemented within this framework. The complete declaration is:

      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-----------------------------------------------------------------------
C     UMAT for single crystal plasticity
C     Developed by Yonggang Huang
C-----------------------------------------------------------------------
      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)
C     ... user coding ...
      RETURN
      END

3.2 变量说明表 / Variable Specification Table 3.2 Variable Specification Table

下表列出了UMAT接口中与单晶塑性实现密切相关的输入/输出变量：

The following table lists the input/output variables in the UMAT interface most relevant to the single-crystal plasticity implementation:

3.3 STATEV布局 / State Variable Layout Table 3.3 State Variable Layout Table

STATEV数组是UMAT子程序的核心数据结构，用于在增量步之间传递和保存材料的历史状态。 黄永刚UMAT的STATEV布局如下（以FCC 12滑移系为例）：

The STATEV array is the core data structure of the UMAT subroutine, used to pass and preserve material history states between increments. The STATEV layout in Huang’s UMAT is as follows (for FCC 12 slip systems as an example):

04 材料参数 / Material Parameters Material Parameters

黄永刚UMAT通过PROPS数组传入全部材料参数。PROPS的索引布局经过精心设计， 将弹性常数、滑移系定义、晶体取向、率相关参数、硬化参数和积分控制参数分区存储。 以下各节详细列出每个PROPS索引对应的参数含义。

Huang’s UMAT receives all material parameters through the PROPS array. The PROPS index layout is carefully designed, partitioning elastic constants, slip system definitions, crystal orientation, rate-dependent parameters, hardening parameters, and integration control parameters. The following sections detail the meaning of each PROPS index.

4.1 弹性常数 PROPS(1-21) / Elastic Constants 4.1 Elastic Constants PROPS(1-21)

弹性常数区定义晶体材料的各向异性弹性行为。对于立方晶系，只需要3个独立常数（C11, C12, C44）， 但程序接口预留了完整21个参数以支持一般各向异性材料（三斜晶系）。

The elastic constant section defines the anisotropic elastic behavior of the crystal material. For cubic crystals, only 3 independent constants (C11, C12, C44) are needed, but the program interface reserves a full 21 parameters to support general anisotropic materials (triclinic).

4.2 滑移系定义 PROPS(25-56) / Slip System Definition 4.2 Slip System Definition PROPS(25-56)

滑移系定义区用于指定每个滑移系的滑移方向s和滑移面法向m（在晶体坐标系中）。 每组滑移系需要3个方向分量和3个法向分量。对于FCC结构，通常定义12个滑移系（{111}<110>）。

The slip system definition section specifies the slip direction s and slip plane normal m (in crystal coordinates) for each slip system. Each slip system requires 3 direction components and 3 normal components. For FCC, 12 slip systems ({111}<110>) are typically defined.

FCC晶体典型的12个滑移系为{111}<110>族，分为4个滑移面，每个面有3个滑移方向。 以下为晶体坐标系中的标准定义（无需归一化，程序内部会自动处理）：

The typical 12 slip systems for FCC crystals belong to the {111}<110> family, divided into 4 slip planes with 3 slip directions each. The following is the standard definition in crystal coordinates (normalization is handled internally by the program):

C     FCC {111}<110> slip systems (standard definition)
C     Plane 1 (111): directions [1-10], [10-1], [01-1]
C     Plane 2 (1-11): directions [110], [101], [0-11]
C     Plane 3 (-111): directions [110], [10-1], [011]
C     Plane 4 (11-1): directions [1-10], [101], [011]

4.3 晶体取向 PROPS(57-70) / Crystal Orientation 4.3 Crystal Orientation PROPS(57-70)

晶体取向区定义晶体坐标系与全局（试样）坐标系之间的旋转关系。 通常采用Bunge Euler角(φ1, Φ, φ2)或方向余弦矩阵表示。

The crystal orientation section defines the rotation relationship between the crystal coordinate system and the global (sample) coordinate system. Bunge Euler angles (φ1, Φ, φ2) or the direction cosine matrix are commonly used.

4.4 率相关参数 PROPS(73-96) / Rate-Dependent Parameters 4.4 Rate-Dependent Parameters PROPS(73-96)

率相关参数区控制幂律流动法则的行为。对于多族滑移系，每族可以独立指定率相关参数。

The rate-dependent parameter section controls the behavior of the power-law flow rule. For multiple families of slip systems, each family can have independent rate-dependent parameters.

4.5 硬化参数 PROPS(97-144) / Hardening Parameters 4.5 Hardening Parameters PROPS(97-144)

硬化参数区包含自硬化和潜硬化的全部参数。根据选择的硬化模型类型（PROPS中的标志），这些参数的解释有所不同。

The hardening parameter section contains all parameters for self-hardening and latent hardening. Their interpretation depends on the selected hardening model type (flag in PROPS).

4.6 积分控制 PROPS(145-160) / Integration Control 4.6 Integration Control PROPS(145-160)

积分控制参数影响隐式Newton-Raphson迭代的数值行为，包括收敛容差、最大迭代次数、 自动时间步长控制策略等。

Integration control parameters affect the numerical behavior of the implicit Newton-Raphson iteration, including convergence tolerance, maximum iteration count, and automatic time step control strategy.

05 代码结构 / Code Structure Code Structure

5.1 子程序调用关系图 / Subroutine Call Graph 5.1 Subroutine Call Graph

黄永刚UMAT由主程序UMAT和多个辅助子程序组成，各子程序承担明确的职责。 以下为调用关系图：

Huang’s UMAT consists of the main UMAT subroutine and multiple auxiliary subroutines, each with a well-defined responsibility. The call graph is as follows:

5.2 主算法流程图 / Main Algorithm Flowchart 5.2 Main Algorithm Flowchart

每个增量步中，UMAT主程序执行以下算法流程：

In each increment, the main UMAT routine executes the following algorithm:

5.3 关键代码段 / Key Code Segments 5.3 Key Code Segments

5.3.1 弹性试应力计算 / Elastic Trial Stress 5.3.1 Elastic Trial Stress

C-----------------------------------------------------------------------
C     弹性试应力计算 / Elastic trial stress computation
C     基于增量开始时的Fe和整个应变增量计算试应力
C-----------------------------------------------------------------------
      DO I=1,3
         DO J=1,3
            FE_INV(I,J)=DFGRD1(I,J)   ! 初始估计: F^e = F (忽略塑性)
         ENDDO
      ENDDO
C
C     从STATEV恢复上一增量步的F^p和F^e
      CALL KINTOTAL(STATEV,NSTATV,FE_FP,NSLIP)
C
C     计算试弹性变形梯度: F^{e,trial} = F_{n+1} * (F^p_n)^{-1}
      CALL LUDCMP(FP_TRY,3,3,INDX,D)
      CALL LUBKSB(FP_TRY,3,3,INDX,FE_TRY_COL1)
      CALL LUBKSB(FP_TRY,3,3,INDX,FE_TRY_COL2)
      CALL LUBKSB(FP_TRY,3,3,INDX,FE_TRY_COL3)
C
C     计算Green弹性应变 E^e = 0.5*(F^{eT}F^e - I)
      DO I=1,3
         DO J=1,3
            FE_TF_FE(I,J)=0.0
            DO K=1,3
               FE_TF_FE(I,J)=FE_TF_FE(I,J)+FE_TRY(K,I)*FE_TRY(K,J)
            ENDDO
         ENDDO
      ENDDO
      DO I=1,3
         DO J=1,3
            EE(I,J)=0.5*(FE_TF_FE(I,J)-DELTA(I,J))
         ENDDO
      ENDDO

5.3.2 流动律实现 / Flow Law Implementation 5.3.2 Flow Law Implementation

C-----------------------------------------------------------------------
C     率相关流动律计算 / Rate-dependent flow rule computation
C     gamma_dot = gamma0 * |tau/g|^n * sign(tau/g)
C-----------------------------------------------------------------------
      SUBROUTINE STRAINRATE(TAU, G, GAMMA0, AM, NSLIP, DGAMMA)
      IMPLICIT NONE
      INTEGER NSLIP, I
      REAL*8 TAU(NSLIP), G(NSLIP), GAMMA0, AM
      REAL*8 DGAMMA(NSLIP), RATIO, ABSRAT, POWER
C
      DO I=1,NSLIP
         IF (G(I).LE.0.0) THEN
            WRITE(*,*) 'ERROR: g(',I,') =',G(I)
            STOP
         ENDIF
         RATIO = TAU(I)/G(I)
         ABSRAT = DABS(RATIO)
C        避免数值溢出 / Avoid numerical overflow
         IF (ABSRAT.GT.1.0D+10) THEN
            POWER = 1.0D+10**AM
         ELSE
            POWER = ABSRAT**AM
         ENDIF
C        剪切应变率 / Shear strain rate
         DGAMMA(I) = GAMMA0 * POWER
         IF (RATIO.LT.0.0) DGAMMA(I) = -DGAMMA(I)
      ENDDO
      RETURN
      END

5.3.3 自硬化实现 / Self-Hardening Implementation 5.3.3 Self-Hardening Implementation

C-----------------------------------------------------------------------
C     双曲正割自硬化模型 / Hyperbolic secant self-hardening model
C     h(gamma) = H0 * sech^2( H0*gamma / (tau_s - tau_0) )
C-----------------------------------------------------------------------
      FUNCTION HSELF(GAMMA, H0, TAUS, TAU0)
      IMPLICIT NONE
      REAL*8 HSELF, GAMMA, H0, TAUS, TAU0
      REAL*8 ARG, SECH2
C
      IF (TAUS.LE.TAU0) THEN
         HSELF = 0.0
         RETURN
      ENDIF
C
      ARG = H0 * GAMMA / (TAUS - TAU0)
C     sech^2(x) = 1/cosh^2(x)
      SECH2 = 1.0 / (DCOSH(ARG)**2)
      HSELF = H0 * SECH2
      RETURN
      END
C
C-----------------------------------------------------------------------
C     当前滑移系强度 / Current slip system strength
C     g(gamma) = tau_s + (tau_0 - tau_s) * tanh( H0*gamma/(tau_s-tau_0) )
C-----------------------------------------------------------------------
      FUNCTION GSLP0(GAMMA, H0, TAUS, TAU0)
      IMPLICIT NONE
      REAL*8 GSLP0, GAMMA, H0, TAUS, TAU0
      REAL*8 ARG
C
      IF (TAUS.LE.TAU0) THEN
         GSLP0 = TAU0
         RETURN
      ENDIF
C
      ARG = H0 * GAMMA / (TAUS - TAU0)
      GSLP0 = TAUS + (TAU0 - TAUS) * DTANH(ARG)
      RETURN
      END

5.3.4 Newton-Raphson迭代核心 / NR Iteration Core 5.3.4 NR Iteration Core

C-----------------------------------------------------------------------
C     Newton-Raphson迭代求解剪切应变增量 / NR iteration for shear increment
C     未知量: DGA(NSLIP) = Δγ^(α)
C     残差:   R(α) = Δγ^(α) - Δt * gamma_dot^{n+1}_(α) = 0
C-----------------------------------------------------------------------
      SUBROUTINE ITERATION(...)
      IMPLICIT NONE
      INTEGER NSLIP, ITMAX, ITER, I, J, INFO
      REAL*8 TOL, DTIME, DGA(NSLIP), GA(NSLIP), TAU(NSLIP)
      REAL*8 G(NSLIP), DG(NSLIP), GAMMA0, AM
      REAL*8 JAC(NSLIP,NSLIP), RES(NSLIP), DGA_INC(NSLIP)
      REAL*8 DFDX(NSLIP), DHDG(NSLIP,NSLIP)
      REAL*8 DET, RESNORM, RESNORM0
      INTEGER INDX(NSLIP)
C
      RESNORM0 = 1.0
      DO ITER=1,ITMAX
C        计算当前应力状态 / Compute current stress state
         CALL STRESSUPDATE(DGA, TAU, SIGMA, ...)
C
C        计算流动律 / Compute flow rule
         CALL STRAINRATE(TAU, G, GAMMA0, AM, NSLIP, DG)
C
C        计算残差 / Compute residual
         DO I=1,NSLIP
            RES(I) = DGA(I) - DTIME * DG(I)
         ENDDO
C
C        计算残差范数 / Compute residual norm
         RESNORM = 0.0
         DO I=1,NSLIP
            RESNORM = RESNORM + RES(I)*RES(I)
         ENDDO
         RESNORM = DSQRT(RESNORM)
         IF (ITER.EQ.1) RESNORM0 = RESNORM
C
C        检查收敛 / Check convergence
         IF (RESNORM.LT.TOL*RESNORM0 .OR. RESNORM.LT.1.0D-12) THEN
            RETURN   ! 收敛 / Converged
         ENDIF
C
C        计算Jacobian矩阵 / Compute Jacobian matrix
C        J_{αβ} = δ_{αβ} - Δt * (∂gamma_dot_α/∂DGA_β)
         DO I=1,NSLIP
            DO J=1,NSLIP
               JAC(I,J) = -DTIME * DHDG(I,J)
            ENDDO
            JAC(I,I) = JAC(I,I) + 1.0
         ENDDO
C
C        LU分解并求解 / LU decomposition and solve
         CALL LUDCMP(JAC, NSLIP, NSLIP, INDX, DET)
         DO I=1,NSLIP
            DGA_INC(I) = -RES(I)
         ENDDO
         CALL LUBKSB(JAC, NSLIP, NSLIP, INDX, DGA_INC)
C
C        更新未知量 / Update unknowns
         DO I=1,NSLIP
            DGA(I) = DGA(I) + DGA_INC(I)
         ENDDO
      ENDDO
C
      WRITE(*,*) 'WARNING: NR iteration not converged in UMAT'
      RETURN
      END

5.3.5 DDSDDE一致切线刚度 / Consistent Tangent DDSDDE 5.3.5 Consistent Tangent DDSDDE

C-----------------------------------------------------------------------
C     计算一致切线刚度矩阵 DDSDDE / Compute consistent tangent modulus
C     DDSDDE = ∂Δσ / ∂Δε
C-----------------------------------------------------------------------
      SUBROUTINE CONSISTENTTANGENT(...)
      IMPLICIT NONE
      INTEGER NTENS, NSLIP, I, J, K, L
      REAL*8 DDSDDE(NTENS,NTENS), CEL(6,6)
      REAL*8 P_ALPHA(3,3,NSLIP), DGA(NSLIP)
      REAL*8 DTAU_DSIG(NTENS,NSLIP), DHDG(NSLIP,NSLIP)
      REAL*8 JAC(NSLIP,NSLIP), INV_JAC(NSLIP,NSLIP)
      REAL*8 TERM(NTENS,NSLIP), SUM
      INTEGER INDX(NSLIP)
      REAL*8 DET
C
C     弹性刚度矩阵 / Elastic stiffness matrix
      DO I=1,NTENS
         DO J=1,NTENS
            DDSDDE(I,J) = CEL(I,J)
         ENDDO
      ENDDO
C
C     计算 ∂tau_α/∂σ 并组装修正项
C     Compute ∂tau_α/∂σ and assemble correction
      DO K=1,NSLIP
         DO I=1,NTENS
            TERM(I,K) = 0.0
            DO J=1,NTENS
               TERM(I,K) = TERM(I,K) + CEL(I,J)*P_ALPHA_VOIGT(J,K)
            ENDDO
         ENDDO
      ENDDO
C
C     计算 (I - Δt * ∂dg/∂g * ∂g/∂DGA)^{-1}
      CALL LUDCMP(JAC, NSLIP, NSLIP, INDX, DET)
      CALL MATINV(JAC, NSLIP, NSLIP, INDX, INV_JAC)
C
C     DDSDDE = C^e - C^e:P * INV_JAC * DTIME * dg/dtau * d tau/d eps
      DO I=1,NTENS
         DO J=1,NTENS
            SUM = 0.0
            DO K=1,NSLIP
               DO L=1,NSLIP
                  SUM = SUM + TERM(I,K)*INV_JAC(K,L)*DTAU_DSIG(L,J)
               ENDDO
            ENDDO
            DDSDDE(I,J) = DDSDDE(I,J) - SUM * DTIME
         ENDDO
      ENDDO
      RETURN
      END

06 Jacobian矩阵 / Jacobian Matrix Jacobian Matrix

6.1 应变分解 / Strain Decomposition 6.1 Strain Decomposition

在Abaqus/Standard的隐式求解框架中，UMAT必须提供一致切线刚度矩阵DDSDDE， 即应力增量对应变增量的偏导数DDSDDE = ∂Δσ / ∂Δε。 对于弹塑性材料，总应变增量分解为弹性部分和塑性部分：

In the Abaqus/Standard implicit solver framework, UMAT must provide the consistent tangent stiffness matrix DDSDDE, i.e., the partial derivative of stress increment with respect to strain increment, DDSDDE = ∂Δσ / ∂Δε. For elastoplastic materials, the total strain increment is decomposed into elastic and plastic parts:

在小变形假设下，弹性应变增量与应力增量的关系由广义胡克定律给出：

Under the small-deformation assumption, the relationship between elastic strain increment and stress increment is given by the generalized Hooke’s law:

塑性应变增量由各滑移系的剪切应变增量贡献：

The plastic strain increment is contributed by the shear strain increments on all slip systems:

因此，一致切线刚度矩阵的关键在于计算 ∂Δγ^(α) / ∂Δε， 即剪切应变增量对应变增量的敏感度。

Therefore, the key to the consistent tangent stiffness matrix is computing ∂Δγ^(α) / ∂Δε, i.e., the sensitivity of shear strain increment to strain increment.

6.2 线性方程组 / Linear System 6.2 Linear System

在Newton-Raphson迭代中，每个滑移系的剪切应变增量满足以下非线性方程：

In the Newton-Raphson iteration, the shear strain increment for each slip system satisfies the following nonlinear equation:

对上述方程关于Δγ^(β)求偏导，得到Newton-Raphson迭代的Jacobian矩阵（局部Newton迭代矩阵）：

Taking the partial derivative of the above equation with respect to Δγ^(β) yields the Jacobian matrix for Newton-Raphson iteration (local Newton iteration matrix):

各项偏导数的物理含义和表达式如下：

The physical meaning and expressions of each partial derivative are:

• ∂γ̇^(α) / ∂τ^(α): 流动律对分解剪应力的敏感度，等于 n · γ̇^(α) / τ^(α) Sensitivity of flow rule to resolved shear stress, equal to n · γ̇^(α) / τ^(α)
• ∂τ^(α) / ∂Δγ^(β): 分解剪应力对剪切增量的敏感度，等于 -C^e : P^(β) : P^(α)（应力松弛效应） Sensitivity of resolved shear stress to shear increment, equal to -C^e : P^(β) : P^(α) (stress relaxation effect)
• ∂g^(α) / ∂Δγ^(β): 滑移系强度对剪切增量的敏感度，等于 h_αβ · sign(τ^(β)) Sensitivity of slip system strength to shear increment, equal to h_αβ · sign(τ^(β))
• ∂γ̇^(α) / ∂g^(α): 流动律对强度的敏感度，等于 -n · γ̇^(α) / g^(α) Sensitivity of flow rule to strength, equal to -n · γ̇^(α) / g^(α)

将所有偏导数代入后，局部Newton矩阵可显式写出。该矩阵的大小为NSLIP × NSLIP（例如FCC为12×12）。 在每次Newton迭代中，对该矩阵进行LU分解并回代求解剪切应变增量修正值ΔDGA。

After substituting all partial derivatives, the local Newton matrix can be written explicitly. Its size is NSLIP × NSLIP (e.g., 12×12 for FCC). In each Newton iteration, this matrix undergoes LU decomposition and back-substitution to solve for the shear strain increment correction ΔDGA.

6.3 DDSDDE构造 / DDSDDE Construction 6.3 DDSDDE Construction

一致切线刚度矩阵的构造公式为：

The consistent tangent stiffness matrix is constructed as:

上述公式中，第一项C^e是纯弹性刚度，第二项是塑性修正。 当材料处于纯弹性状态（无滑移激活）时，第二项为零，DDSDDE退化为弹性刚度矩阵。 当大量滑移系激活时，第二项显著降低材料的表观刚度，反映塑性流动引起的应力松弛。

In the above formula, the first term C^e is the pure elastic stiffness, and the second term is the plastic correction. When the material is in a purely elastic state (no slip activated), the second term vanishes and DDSDDE reduces to the elastic stiffness matrix. When many slip systems are active, the second term significantly reduces the apparent material stiffness, reflecting stress relaxation due to plastic flow.

具体计算步骤如下：

The specific calculation steps are:

1. 计算弹性刚度矩阵C^e的Voigt形式（6×6矩阵） Compute the Voigt form of elastic stiffness matrix C^e (6×6 matrix)
2. 计算每个滑移系的Schmid张量P^(α)的Voigt向量形式 Compute the Voigt vector form of Schmid tensor P^(α) for each slip system
3. 计算流动律导数 ∂γ̇^(α) / ∂τ^(α) = n · γ̇^(α) / τ^(α) Compute flow rule derivative ∂γ̇^(α) / ∂τ^(α) = n · γ̇^(α) / τ^(α)
4. 组装局部Newton矩阵J_αβ并进行LU分解求逆 Assemble local Newton matrix J_αβ and perform LU decomposition for inversion
5. 计算中间张量：对每个滑移系α，计算 Q^(α)_ij = (J^-1_αβ · DFLOW_β) · (P^(β)_kl C^e_klmn) Compute intermediate tensor: for each slip system α, compute Q^(α)_ij = (J^-1_αβ · DFLOW_β) · (P^(β)_kl C^e_klmn)
6. 组装DDSDDE：DDSDDE_ij = C^e_ij – Σ_α (C^e_ijkl P^(α)_kl) Q^(α)_mn Assemble DDSDDE: DDSDDE_ij = C^e_ij – Σ_α (C^e_ijkl P^(α)_kl) Q^(α)_mn

DDSDDE的对称性取决于材料模型。对于关联流动法则和二次屈服面，DDSDDE是对称的； 但对于非关联流动或复杂的潜硬化模型，DDSDDE可能非对称。 黄永刚UMAT中，由于采用关联的Schmid流动和一般的硬化矩阵，DDSDDE通常具有轻微的非对称性。 Abaqus/Standard支持非对称刚度矩阵，但需要在*Step中设置UNSymm=YES以优化求解器性能。

The symmetry of DDSDDE depends on the material model. For associative flow rules and quadratic yield surfaces, DDSDDE is symmetric; but for non-associative flow or complex latent hardening models, DDSDDE may be non-symmetric. In Huang’s UMAT, due to the associative Schmid flow and general hardening matrix, DDSDDE usually has slight non-symmetry. Abaqus/Standard supports non-symmetric stiffness matrices, but *Step requires UNSymm=YES to optimize solver performance.

07 典型参数 / Typical Parameters Typical Parameters

以下给出FCC铜和FCC铝的典型材料参数，这些参数在大量文献和实验中得到验证， 可作为Abaqus输入文件中*User Material的参考值。

The following presents typical material parameters for FCC copper and FCC aluminum, validated in numerous literature and experiments, serving as reference values for the *User Material keyword in Abaqus input files.

7.1 FCC铜 / FCC Copper 7.1 FCC Copper

铜是FCC结构，室温下主要滑移系为{111}<110>，共12个。 以下参数适用于室温（约300K）下的单晶铜模拟，采用双曲正割硬化模型。

Copper has an FCC structure. At room temperature, the primary slip systems are {111}<110>, totaling 12. The following parameters are suitable for single-crystal copper simulation at room temperature (~300K) using the hyperbolic secant hardening model.

7.2 FCC铝 / FCC Aluminum 7.2 FCC Aluminum

铝同样是FCC结构，但其层错能较高，位错交滑移更容易发生，因此硬化行为和铜有明显差异。 以下参数适用于高纯铝单晶在室温下的模拟。

Aluminum is also FCC, but with higher stacking fault energy, making dislocation cross-slip easier, so its hardening behavior differs significantly from copper. The following parameters are suitable for high-purity single-crystal aluminum at room temperature.

08 使用注意事项 / Usage Notes Usage Notes

以下列出使用黄永刚UMAT子程序时常见的注意事项和最佳实践建议， 涵盖输入文件设置、材料参数校准、收敛性调试和后处理等方面。

The following lists common considerations and best-practice recommendations when using Huang’s UMAT subroutine, covering input file setup, material parameter calibration, convergence debugging, and post-processing.

8.1 Abaqus输入文件设置 / Input File Setup 8.1 Input File Setup

• *User Material关键字：必须在输入文件中通过*User Material, constants=NPROPS定义PROPS数组长度，并确保*Depvar的NSTATV足够大以容纳所有状态变量。对于FCC 12滑移系，建议NSTATV >= 70。 *User Material keyword: The PROPS array length must be defined via *User Material, constants=NPROPS, and *Depvar’s NSTATV must be large enough for all state variables. For FCC 12 slip systems, NSTATV >= 70 is recommended.
• UNSymm选项：由于晶体塑性的一致切线刚度矩阵通常轻微非对称，建议在*Step中使用UNSymm=YES以启用非对称求解器，提高收敛速度和稳定性。 UNSymm option: Because the consistent tangent stiffness matrix of crystal plasticity is usually slightly non-symmetric, it is recommended to use UNSymm=YES in *Step to enable the non-symmetric solver, improving convergence speed and stability.
• 时间步长控制：晶体塑性UMAT对时间步长敏感。建议在初始加载阶段使用较小的时间步，待塑性变形稳定后可适度增大。可通过*Static中的INITIAL和MAXIMUM参数控制。 Time step control: Crystal plasticity UMAT is sensitive to time step size. It is recommended to use small time steps in the initial loading stage, and moderately increase after plastic deformation stabilizes. Control via INITIAL and MAXIMUM in *Static.
• 单元类型选择：推荐使用C3D8R（8节点六面体减缩积分单元）或C3D8（完全积分）。避免使用线性三角形/四面体单元模拟大应变局部化，因为它们在不可压缩塑性变形中会产生体积锁定。 Element type selection: C3D8R (8-node hexahedron with reduced integration) or C3D8 (full integration) are recommended. Avoid linear triangular/tetrahedral elements for large-strain localization, as they suffer volumetric locking in incompressible plastic deformation.
• 沙漏控制：使用C3D8R单元时，必须配合*Hourglass Stiffness或*Section Controls设置合理的沙漏刚度，以抑制零能模式。 Hourglass control: When using C3D8R elements, *Hourglass Stiffness or *Section Controls must be set with reasonable hourglass stiffness to suppress zero-energy modes.

8.2 材料参数校准 / Parameter Calibration 8.2 Parameter Calibration

• 弹性常数：C11、C12、C44应从超声共振或单晶弹性实验获取。注意Abaqus中应力单位为MPa，弹性常数也必须以MPa输入。 Elastic constants: C11, C12, C44 should be obtained from ultrasonic resonance or single-crystal elastic experiments. Note that Abaqus uses MPa for stress, so elastic constants must also be entered in MPa.
• 初始CRSS tau0：可通过单晶初始屈服应力结合Schmid因子估算。对于不同取向的单晶试样，tau0应唯一，而宏观屈服应力随取向变化。 Initial CRSS tau0: Can be estimated from single-crystal initial yield stress combined with the Schmid factor. For single-crystal specimens of different orientations, tau0 should be unique, while macroscopic yield stress varies with orientation.
• 率敏感指数n：应通过不同应变速率下的单晶拉伸实验确定。n值对流动应力水平有显著影响，但对硬化模量影响较小。若实验数据不足，可参考文献中同种材料的推荐值。 Rate sensitivity exponent n: Should be determined from single-crystal tensile experiments at different strain rates. n significantly affects flow stress level but has minor effect on hardening modulus. If experimental data is insufficient, refer to recommended values in the literature for the same material.
• 硬化参数：H0、taus和tau0应通过单晶应力-应变曲线的塑性段拟合确定。建议使用非线性最小二乘法对实验数据进行优化拟合。 Hardening parameters: H0, taus, and tau0 should be determined by fitting the plastic portion of single-crystal stress-strain curves. Nonlinear least-squares optimization is recommended for fitting experimental data.
• 潜硬化系数q：各向同性近似下可取q=1.0。若要更精确描述，需要通过双滑移或多滑移实验标定Bassani模型中的7个潜硬化系数。 Latent hardening coefficient q: In isotropic approximation, q=1.0 can be used. For more accurate description, the 7 latent hardening coefficients in the Bassani model need to be calibrated through double-slip or multi-slip experiments.
• 晶体取向：单晶模拟中取向参数至关重要。建议通过EBSD实验测量实际试样的Euler角，而非仅依赖名义取向。即使微小的取向偏差也可能导致宏观应力-应变响应的显著差异。 Crystal orientation: Orientation parameters are crucial in single-crystal simulations. It is recommended to measure the Euler angles of actual specimens via EBSD, rather than relying solely on nominal orientations. Even small orientation deviations can cause significant differences in macroscopic stress-strain response.

8.3 收敛性调试 / Convergence Debugging 8.3 Convergence Debugging

• 时间步过大：若出现”Too many attempts made for this increment”错误，最常见的原因是时间步过大导致局部Newton迭代不收敛。可通过减小最大时间步或启用自动步长控制解决。 Time step too large: If “Too many attempts made for this increment” appears, the most common cause is excessive time step size causing local Newton iteration divergence. Reduce maximum time step or enable automatic step control.
• 奇异Jacobian：当多个滑移系同时激活且硬化矩阵接近奇异时，LU分解可能失败。可通过增加收敛容差TOL或在流动律中引入正则化来改善。 Singular Jacobian: When multiple slip systems activate simultaneously and the hardening matrix is nearly singular, LU decomposition may fail. Increase convergence tolerance TOL or introduce regularization in the flow rule.
• 负剪切率：若某些滑移系的分解剪应力为负且绝对值很大，流动律可能产生不合理的负剪切率。程序中应有符号处理逻辑，确保物理合理性。 Negative shear rate: If resolved shear stress on some slip systems is negative with large magnitude, the flow rule may produce unreasonable negative shear rates. The program should have sign-handling logic to ensure physical consistency.
• 检查STATEV初始化：STATEV的初始值必须正确设置，特别是弹性变形梯度Fe应初始化为单位矩阵，Fp也初始化为单位矩阵，g初始化为tau0。任何错误的初值都可能导致第一步即发散。 Check STATEV initialization: Initial STATEV values must be set correctly, especially Fe and Fp initialized to identity matrices, and g initialized to tau0. Any incorrect initial value may cause divergence in the first step.
• PNEWDT的使用：当局部迭代超过ITMAX仍未收敛时，应在UMAT中设置PNEWDT < 1（如0.5），告知Abaqus减小当前增量步并重试。这是一种有效的自适应步长策略。 Use of PNEWDT: When local iteration exceeds ITMAX without convergence, set PNEWDT < 1 (e.g., 0.5) in UMAT to instruct Abaqus to reduce the current increment and retry. This is an effective adaptive step strategy.
• 输出调试信息：设置OUTPUT=2可启用详细调试输出，包括每次Newton迭代的残差范数、剪切增量和应力状态。这对于定位发散原因非常有帮助。 Output debug information: Setting OUTPUT=2 enables verbose debug output, including residual norm, shear increments, and stress state for each Newton iteration. This is very helpful for locating divergence causes.

8.4 后处理与结果验证 / Post-Processing and Verification 8.4 Post-Processing and Verification

• 应力-应变曲线：将Abaqus输出的宏观应力-应变曲线与单晶实验数据进行对比，验证模型预测的基本力学行为是否正确。 Stress-strain curve: Compare the macroscopic stress-strain curve from Abaqus with single-crystal experimental data to verify whether the basic mechanical behavior predicted by the model is correct.
• 滑移系激活分析：通过后处理STATEV中存储的累积剪切应变，可以识别哪些滑移系被激活及其相对活动强度。这对于理解变形机制和织构演化至关重要。 Slip system activation analysis: By post-processing accumulated shear strains stored in STATEV, one can identify which slip systems are activated and their relative activity levels. This is crucial for understanding deformation mechanisms and texture evolution.
• 织构验证：如果启用了取向更新（PROPS(70)=1），应提取变形后的晶体取向并与实验EBSD织构数据进行对比，评估织构演化预测的准确性。 Texture verification: If orientation update is enabled (PROPS(70)=1), extract deformed crystal orientations and compare with experimental EBSD texture data to assess the accuracy of texture evolution prediction.
• 网格收敛性：进行网格敏感性分析，确保结果不依赖于离散化程度。晶体塑性模拟中，应变局部化带（shear band）的宽度通常受网格尺寸影响，应使用足够细密的网格。 Mesh convergence: Perform mesh sensitivity analysis to ensure results do not depend on discretization level. In crystal plasticity simulations, the width of strain localization bands (shear bands) is usually affected by mesh size; sufficiently fine meshes should be used.
• 能量守恒检查：验证塑性耗散功与外部功的关系，确保数值实现没有能量错误。Abaqus输出的能量历史（ALLIE、ALLKE、ALLPD等）可用于此目的。 Energy conservation check: Verify the relationship between plastic dissipation work and external work to ensure the numerical implementation has no energy errors. Abaqus energy history outputs (ALLIE, ALLKE, ALLPD, etc.) can be used for this purpose.

8.5 常见错误与解决方案 / Common Errors and Solutions 8.5 Common Errors and Solutions

Ref 参考文献 / References References

以下列出黄永刚UMAT子程序相关的核心理论文献，按发表时间排序：

The following lists core theoretical references related to Huang’s UMAT subroutine, in chronological order:

App 符号表 / Notation Notation

本文档中使用的关键数学符号及其含义如下表所示：

The key mathematical symbols used in this document and their meanings are listed in the following table: