Chapter 1 概述Chapter 1 Overview

1.1 VUMAT简介1.1 Introduction to VUMAT

VUMAT（Vectorized User Material）是Abaqus/Explicit求解器中用于自定义材料本构的子程序接口。与适用于隐式求解器Abaqus/Standard的UMAT不同，VUMAT针对显式动力学分析进行了专门优化，采用批量处理（NBLOCK）模式一次性处理多个积分点，从而充分利用向量化和并行计算能力。VUMAT (Vectorized User Material) is the subroutine interface in Abaqus/Explicit for defining custom constitutive material behavior. Unlike UMAT, which is designed for the implicit solver Abaqus/Standard, VUMAT is specifically optimized for explicit dynamic analysis. It processes multiple integration points in batch mode (NBLOCK), thereby fully utilizing vectorization and parallel computing capabilities.

黄永刚教授开发的晶体塑性VUMAT子程序，将单晶塑性理论拓展至显式有限元框架。该子程序基于率相关晶体塑性本构模型，通过显式时间积分实现高效的大规模晶体塑性有限元模拟，特别适用于冲击、爆炸、高速变形等强非线性问题。Professor Yonggang Huang’s crystal plasticity VUMAT subroutine extends single-crystal plasticity theory into the explicit finite element framework. Based on a rate-dependent crystal plasticity constitutive model, it achieves efficient large-scale crystal plasticity finite element simulations through explicit time integration, making it particularly suitable for highly nonlinear problems such as impact, explosion, and high-speed deformation.

晶体塑性有限元方法（CPFEM）将宏观有限元与微观晶体学相结合，通过追踪每个积分点上数千个晶粒的取向演化，揭示材料变形过程中的织构演化、应力集中和局部化现象。VUMAT版本保留了UMAT版本的核心物理模型，但在计算架构上进行了根本性重构。Crystal Plasticity Finite Element Method (CPFEM) combines macroscopic finite elements with microscopic crystallography. By tracking the orientation evolution of thousands of grains at each integration point, it reveals texture evolution, stress concentration, and localization phenomena during material deformation. The VUMAT version retains the core physical model of the UMAT version but undergoes fundamental restructuring in computational architecture.

1.2 与UMAT的关系1.2 Relationship with UMAT

黄永刚VUMAT与UMAT共享相同的物理本构核心：两者都基于Schmid定律描述滑移系上的分切应力与剪切应变率关系，都采用Newton-Raphson迭代求解非线性方程组，都使用相同的材料参数定义方式。Huang’s VUMAT and UMAT share the same physical constitutive core: both are based on Schmid’s law to describe the relationship between resolved shear stress and shear strain rate on slip systems; both employ Newton-Raphson iteration to solve nonlinear equation systems; and both use identical material parameter definitions.

主要差异在于接口层和积分层：UMAT针对Standard隐式求解器设计，每个增量步调用一次单个积分点，需要用户提供材料雅可比矩阵DDSDDE；VUMAT针对Explicit显式求解器设计，每个增量步批量调用NBLOCK个积分点，不需要提供雅可比矩阵，但需要用户自行处理旋转张量DROT和状态变量的新旧值传递。The main differences lie in the interface and integration layers: UMAT is designed for the Standard implicit solver, calling a single integration point once per increment, and requires the user to provide the material Jacobian matrix DDSDDE. VUMAT is designed for the Explicit solver, calling NBLOCK integration points in batch per increment. It does not require the Jacobian matrix but requires the user to handle the rotation tensor DROT and the transfer of old/new state variable values.

1.3 适用场景1.3 Application Scenarios

• 高速冲击与碰撞模拟High-speed impact and collision simulations
• 爆炸与爆轰加载分析Explosion and detonation loading analysis
• 金属切削与成形加工Metal cutting and forming processes
• 动态织构演化研究Dynamic texture evolution studies
• 多晶材料大变形行为Large deformation behavior of polycrystalline materials
• 绝热剪切带形成Adiabatic shear band formation
• 微观组织-宏观性能关联Microstructure-macroscopic property correlations
• 晶体取向对断裂的影响Effect of crystal orientation on fracture

Chapter 2 VUMAT与UMAT对比Chapter 2 VUMAT vs UMAT Comparison

理解VUMAT与UMAT的差异是成功使用黄永刚晶体塑性子程序的关键。以下对比表格详细列出了两个接口在架构、数据传递、计算流程等方面的主要区别。Understanding the differences between VUMAT and UMAT is key to successfully using the Huang crystal plasticity subroutine. The following comparison table details the main differences between the two interfaces in architecture, data transfer, and computational flow.

Chapter 3 VUMAT接口Chapter 3 VUMAT Interface

3.1 完整Fortran签名3.1 Full Fortran Signature

VUMAT子程序的Fortran接口定义如下。与UMAT相比，VUMAT的变量列表更长，反映了显式求解器批量处理的需求。The Fortran interface definition of the VUMAT subroutine is shown below. Compared to UMAT, the VUMAT variable list is longer, reflecting the batch processing requirements of the explicit solver.

      SUBROUTINE VUMAT(
     &  NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, LANNEAL,
     &  STEPTIME, TOTALTIME, DT, CMNAME, COORDMP, CHARLENGTH,
     &  PROPS, DENSITY, STRAININC, RELSPININC,
     &  TEMPOLD, STRETCHOLD, DEFGRADOLD, FIELDOLD,
     &  STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,
     &  TEMPNEW,  STRETCHNEW,  DEFGRADNEW,  FIELDNEW,
     &  STRESSNEW,  STATENEW,  ENERINTERNNEW, ENERINELASNEW )

      CHARACTER*80 CMNAME
      DIMENSION PROPS(NPROPS), DENSITY(NBLOCK), COORDMP(NBLOCK),
     &  STRAININC(NBLOCK,NDIR+NSHR), RELSPININC(NBLOCK,NSHR),
     &  TEMPOLD(NBLOCK), STRETCHOLD(NBLOCK,NDIR+NSHR),
     &  DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR),
     &  FIELDOLD(NBLOCK,NFIELDV), STRESSOLD(NBLOCK,NDIR+NSHR),
     &  STATEOLD(NBLOCK,NSTATEV), ENERINTERNOLD(NBLOCK),
     &  ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),
     &  STRETCHNEW(NBLOCK,NDIR+NSHR),
     &  DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),
     &  FIELDNEW(NBLOCK,NFIELDV), STRESSNEW(NBLOCK,NDIR+NSHR),
     &  STATENEW(NBLOCK,NSTATEV), ENERINTERNNEW(NBLOCK),
     &  ENERINELASNEW(NBLOCK)

      ! User coding to define STRESSNEW and STATENEW
      RETURN
      END

3.2 输入变量表3.2 Input Variables Table

3.3 输出变量表3.3 Output Variables Table

3.4 STATEV布局3.4 State Variable Layout

在黄永刚VUMAT实现中，状态变量数组STATEV按照以下布局存储晶体塑性相关的内部变量。正确理解此布局对于调试和结果提取至关重要。In Huang’s VUMAT implementation, the state variable array STATEV stores crystal-plasticity-related internal variables according to the following layout. Understanding this layout is crucial for debugging and result extraction.

! STATEV layout for Huang VUMAT crystal plasticity
! For FCC crystals with 12 slip systems
!
!  1 -  6:  Cauchy stress components (11,22,33,12,13,23)
!  7 -  9:  Euler angles (phi1, Phi, phi2) in degrees
! 10 - 21:  Shear strain on each slip system (gamma, 12 systems)
! 22 - 33:  Shear strain rate on each slip system (gammadot, 12 systems)
! 34 - 45:  Resolved shear stress on each slip system (tau, 12 systems)
! 46 - 57:  Slip resistance on each slip system (g, 12 systems)
! 58 - 60:  Total strain components (elastic+plastic)
! 61 - 63:  Plastic strain components
! 64:      Accumulated shear strain (sum of |gamma|)
! 65:      Number of Newton iterations (for convergence check)
! 66:      CPU time for constitutive update
! 67+:     Additional user-defined variables

Chapter 4 VUMAT外壳子程序Chapter 4 VUMAT Wrapper

VUMAT外壳子程序（Wrapper）是连接Abaqus/Explicit求解器与核心本构算法的桥梁。它负责处理批量数据格式转换、变量映射、旋转张量计算等任务，然后将每个积分点的数据传递给核心本构子程序。The VUMAT wrapper subroutine is the bridge connecting the Abaqus/Explicit solver with the core constitutive algorithm. It handles batch data format conversion, variable mapping, rotation tensor computation, and then passes each integration point data to the core constitutive subroutine.

4.1 初始化4.1 Initialization

在子程序入口处，首先进行必要的初始化工作，包括声明局部数组、定义滑移系数目、读取材料参数等。At the subroutine entry, necessary initialization is performed, including declaring local arrays, defining the number of slip systems, and reading material parameters.

      SUBROUTINE VUMAT(...)
      ! ... argument list ...
      
      ! Crystal type parameters
      PARAMETER(NSLIP_MAX=12, NSYS_MAX=12)
      
      ! Local arrays for single-point processing
      DIMENSION STRESS(6), STRAN(6), DSTRAN(6)
      DIMENSION STATEV(NSTATEV)
      DIMENSION FOLD(3,3), FNEW(3,3), DROT(3,3)
      DIMENSION PROPS_LOCAL(NPROPS)
      
      ! Copy material properties
      DO K=1,NPROPS
        PROPS_LOCAL(K)=PROPS(K)
      END DO
      
      ! Get crystal type from PROPS(1)
      ! 1=FCC, 2=BCC, 3=HCP
      CRYSTAL_TYPE = NINT(PROPS(1))

4.2 NBLOCK循环4.2 NBLOCK Loop

VUMAT的核心结构是一个遍历所有积分点的循环。在循环内部，对每个积分点执行以下操作：提取旧状态变量、计算旋转张量、转换应变格式、调用核心本构算法、存储新状态变量。The core structure of VUMAT is a loop over all integration points. Inside the loop, for each integration point the following operations are performed: extract old state variables, compute rotation tensor, convert strain format, call the core constitutive algorithm, and store new state variables.

      DO KM=1,NBLOCK
        
        ! Extract old stress
        DO K1=1,NDIR+NSHR
          STRESS(K1)=STRESSOLD(KM,K1)
        END DO
        
        ! Extract old state variables
        DO K1=1,NSTATEV
          STATEV(K1)=STATEOLD(KM,K1)
        END DO
        
        ! Extract deformation gradient
        CALL ExtractDefGrad(DEFGRADOLD,KM,FOLD)
        CALL ExtractDefGrad(DEFGRADNEW,KM,FNEW)
        
        ! Compute rotation increment DROT
        CALL ComputeDROT(FOLD,FNEW,DT,DROT)
        
        ! Convert strain increment
        DO K1=1,NDIR+NSHR
          DSTRAN(K1)=STRAININC(KM,K1)
        END DO
        ! Shear strain conversion for components 4-6
        DO K1=NDIR+1,NDIR+NSHR
          DSTRAN(K1)=DSTRAN(K1)*0.5
        END DO
        
        ! Call core constitutive subroutine
        CALL CRYSTALPLASTICITY(
     &    STRESS,DSTRAN,STATEV,NSTATEV,PROPS_LOCAL,NPROPS,
     &    DROT,NDIR,NSHR,DT,TEMPOLD(KM),TEMPNEW(KM))
        
        ! Store new stress
        DO K1=1,NDIR+NSHR
          STRESSNEW(KM,K1)=STRESS(K1)
        END DO
        
        ! Store new state variables
        DO K1=1,NSTATEV
          STATENEW(KM,K1)=STATEV(K1)
        END DO
        
      END DO

flowchart TD
    A[Start VUMAT
开始VUMAT] --> B[Initialize arrays
初始化数组]
    B --> C{KM = 1 to NBLOCK
循环NBLOCK}
    C --> D[Extract STRESSOLD
提取旧应力]
    D --> E[Extract STATEOLD
提取旧状态]
    E --> F[Extract DEFGRADOLD/NEW
提取变形梯度]
    F --> G[Compute DROT
计算旋转张量]
    G --> H[Convert strain format
转换应变格式]
    H --> I[Call CRYSTALPLASTICITY
调用核心本构]
    I --> J[Store STRESSNEW
存储新应力]
    J --> K[Store STATENEW
存储新状态]
    K --> C
    C -->|Loop end
循环结束| L[Return to Abaqus
返回Abaqus]

4.3 剪切应变转换4.3 Shear Strain Conversion

这是VUMAT与UMAT接口之间最容易出错的地方之一。Abaqus/Explicit传递给VUMAT的剪切应变分量是张量剪切应变的两倍（即工程剪切应变），而核心本构算法通常使用张量剪切应变。因此必须进行转换。This is one of the most error-prone aspects between the VUMAT and UMAT interfaces. The shear strain components passed by Abaqus/Explicit to VUMAT are twice the tensor shear strain (i.e., engineering shear strain), while the core constitutive algorithm typically uses tensor shear strain. Therefore, conversion is necessary.

      ! Abaqus/Explicit provides engineering shear strains
      ! for components 4, 5, 6 in STRAININC
      ! Core algorithm expects tensor shear strains
      ! 
      ! 工程剪切应变 -> 张量剪切应变
      ! Engineering -> Tensor shear strain
      
      DO K1=NDIR+1,NDIR+NSHR
        DSTRAN(K1)=STRAININC(KM,K1)*0.5
      END DO

如果不进行此转换，剪切应力将被错误计算，导致整体应力状态偏离正确解。在调试时，可以通过比较简单剪切变形的解析解来验证转换是否正确。Without this conversion, shear stress will be incorrectly computed, causing the overall stress state to deviate from the correct solution. During debugging, the conversion can be verified by comparing with the analytical solution for simple shear deformation.

4.4 分量顺序交换4.4 Component Order Swap

Abaqus/Explicit中应变和应力的分量顺序为：11, 22, 33, 12, 13, 23。而某些核心本构算法（特别是从UMAT移植过来的代码）可能使用不同的顺序：11, 22, 33, 12, 23, 13。黄永刚原始UMAT代码中使用的是Abaqus标准顺序，但在某些修改版本中需要注意此问题。The component order for strain and stress in Abaqus/Explicit is: 11, 22, 33, 12, 13, 23. Some core constitutive algorithms (especially code ported from UMAT) may use a different order: 11, 22, 33, 12, 23, 13. The original Huang UMAT code uses the Abaqus standard order, but modified versions require attention to this issue.

如果核心子程序期望的顺序与VUMAT传入的不同，需要在调用前后进行分量交换。例如：If the core subroutine expects a different order from what VUMAT passes, component swapping is needed before and after the call. For example:

      ! If core subroutine expects order: 11,22,33,12,23,13
      ! but VUMAT provides: 11,22,33,12,13,23
      ! Swap components 5 and 6 before calling core
      
      TEMP_SWAP=DSTRAN(5)
      DSTRAN(5)=DSTRAN(6)
      DSTRAN(6)=TEMP_SWAP
      
      ! Call core constitutive with swapped order
      CALL CRYSTALPLASTICITY(...)
      
      ! Swap back for output
      TEMP_SWAP=STRESS(5)
      STRESS(5)=STRESS(6)
      STRESS(6)=TEMP_SWAP

4.5 变形梯度展开4.5 Deformation Gradient Expansion

Abaqus/Explicit以压缩数组形式传递变形梯度。对于三维问题，DEFGRAD数组包含9个分量（3×3矩阵），但存储为一维数组。需要正确展开到3×3矩阵以进行后续的张量运算。Abaqus/Explicit passes the deformation gradient in a compressed array form. For 3D problems, the DEFGRAD array contains 9 components (3×3 matrix) but is stored as a one-dimensional array. It must be correctly expanded into a 3×3 matrix for subsequent tensor operations.

      SUBROUTINE ExtractDefGrad(DEFGRAD,KM,F)
      DIMENSION DEFGRAD(NBLOCK,*), F(3,3)
      
      ! Abaqus stores deformation gradient as:
      ! F11, F21, F31, F12, F22, F32, F13, F23, F33
      ! (column-major order for each integration point)
      
      F(1,1)=DEFGRAD(KM,1)
      F(2,1)=DEFGRAD(KM,2)
      F(3,1)=DEFGRAD(KM,3)
      F(1,2)=DEFGRAD(KM,4)
      F(2,2)=DEFGRAD(KM,5)
      F(3,2)=DEFGRAD(KM,6)
      F(1,3)=DEFGRAD(KM,7)
      F(2,3)=DEFGRAD(KM,8)
      F(3,3)=DEFGRAD(KM,9)
      
      RETURN
      END

4.6 变量映射4.6 Variable Mapping

VUMAT外壳的核心任务是将Abaqus/Explicit的批量数组映射到单个积分点的局部变量，并正确处理新旧状态的时间推进。以下表格总结了主要变量的映射关系。The core task of the VUMAT wrapper is to map the batch arrays from Abaqus/Explicit to local variables for a single integration point and correctly handle the time advancement of old/new states. The following table summarizes the mapping of main variables.

Chapter 5 DROT旋转张量计算Chapter 5 DROT Computation

5.1 为什么需要自行计算5.1 Why Self-Computed

在UMAT接口中，Abaqus/Standard会直接将增量旋转张量DROT传递给子程序，用户可以直接使用它来实现应力或取向的客观积分。然而，VUMAT接口不提供DROT变量。这是因为显式求解器的架构不同：它通过变形梯度来描述运动学，要求用户自行从变形梯度中提取旋转信息。In the UMAT interface, Abaqus/Standard directly passes the incremental rotation tensor DROT to the subroutine, which users can directly use for objective integration of stress or orientation. However, the VUMAT interface does not provide the DROT variable. This is because the explicit solver has a different architecture: it describes kinematics through the deformation gradient and requires users to extract rotation information from the deformation gradient themselves.

旋转张量DROT的本质是描述材料点从一个构型到下一个构型的刚性旋转。在晶体塑性中，准确的旋转计算对于追踪晶体取向演化（织构）至关重要。如果DROT计算不准确，会导致取向漂移和应力计算错误。The essence of the rotation tensor DROT is to describe the rigid rotation of a material point from one configuration to the next. In crystal plasticity, accurate rotation computation is crucial for tracking crystal orientation evolution (texture). If DROT is computed inaccurately, it leads to orientation drift and stress calculation errors.

5.2 速度梯度分解5.2 Velocity Gradient from F

从变形梯度计算旋转张量的标准方法是首先计算速度梯度L，然后将其分解为对称部分（变形率D）和反对称部分（旋率W）。The standard method for computing the rotation tensor from the deformation gradient is to first calculate the velocity gradient L, then decompose it into symmetric part (stretching rate D) and antisymmetric part (spin W).

在显式积分中，我们通常使用增量形式。给定增量步开始时的变形梯度F_n和结束时的F_n+1，可以近似计算增量速度梯度。In explicit integration, we typically use incremental form. Given the deformation gradient at the beginning of the increment F_n and at the end F_{{n+1}}, the incremental velocity gradient can be approximately computed.

5.3 旋转增量提取5.3 Rotation Increment Extraction

从旋率W提取旋转增量有多种方法。最简单的是一阶近似：R = I + W·Δt，但这只适用于小旋转。黄永刚VUMAT采用更精确的二阶近似公式，在精度和计算效率之间取得平衡。There are multiple methods to extract the rotation increment from the spin W. The simplest is first-order approximation: R = I + W·Δt, but this is only suitable for small rotations. The Huang VUMAT uses a more accurate second-order approximation formula, balancing accuracy and computational efficiency.

5.4 二阶精度公式5.4 2nd-order Formula

黄永刚VUMAT中使用的DROT计算采用基于旋率张量的二阶精度Cayley变换。对于给定的旋率增量W·Δt，旋转增量张量可以通过以下公式计算：Huang’s VUMAT uses a second-order accurate Cayley transform based on the spin tensor to compute DROT. For a given spin increment W·Δt, the rotation increment tensor can be computed by the following formula:

这个公式等价于对反对称张量W·Δt应用Cayley变换，它保证了旋转张量的正交性（R^TR = I）且行列式为+1，同时具有二阶时间精度。This formula is equivalent to applying the Cayley transform to the antisymmetric tensor W·Δt. It guarantees the orthogonality of the rotation tensor (R^TR = I) with determinant +1, while possessing second-order temporal accuracy.

在实际代码中，为了避免矩阵求逆的数值开销，通常将公式展开为显式分量形式。对于反对称矩阵W，其分量形式为：In actual code, to avoid the numerical cost of matrix inversion, the formula is typically expanded into explicit component form. For the antisymmetric matrix W, its component form is:

      SUBROUTINE ComputeDROT(FOLD,FNEW,DT,DROT)
      IMPLICIT NONE
      DOUBLE PRECISION FOLD(3,3),FNEW(3,3),DT,DROT(3,3)
      DOUBLE PRECISION FINV(3,3),L(3,3),W(3,3)
      DOUBLE PRECISION A(3,3),B(3,3),DET
      DOUBLE PRECISION W1,W2,W3,FAC
      INTEGER I,J,K
      
      ! Compute inverse of FOLD
      CALL M3INV(FOLD,FINV,DET)
      
      ! Compute velocity gradient: L = (FNEW-FOLD)/DT * FINV
      DO I=1,3
        DO J=1,3
          L(I,J)=0.0
          DO K=1,3
            L(I,J)=L(I,J)+(FNEW(I,K)-FOLD(I,K))/DT*FINV(K,J)
          END DO
        END DO
      END DO
      
      ! Extract spin tensor W = 0.5*(L - L^T)
      W1=0.5*(L(3,2)-L(2,3))
      W2=0.5*(L(1,3)-L(3,1))
      W3=0.5*(L(2,1)-L(1,2))
      
      W(1,1)=0.0;   W(1,2)=-W3;   W(1,3)=W2
      W(2,1)=W3;    W(2,2)=0.0;   W(2,3)=-W1
      W(3,1)=-W2;   W(3,2)=W1;    W(3,3)=0.0
      
      ! Second-order approximation: R = (I-0.5*W*dt)^-1 * (I+0.5*W*dt)
      FAC=0.5*DT
      DO I=1,3
        DO J=1,3
          A(I,J)=-FAC*W(I,J)  ! A = -0.5*W*dt
          B(I,J)= FAC*W(I,J)  ! B =  0.5*W*dt
        END DO
        A(I,I)=A(I,I)+1.0     ! A = I - 0.5*W*dt
        B(I,I)=B(I,I)+1.0     ! B = I + 0.5*W*dt
      END DO
      
      CALL M3INV(A,FINV,DET)  ! FINV = A^-1
      
      ! DROT = A^-1 * B
      DO I=1,3
        DO J=1,3
          DROT(I,J)=0.0
          DO K=1,3
            DROT(I,J)=DROT(I,J)+FINV(I,K)*B(K,J)
          END DO
        END DO
      END DO
      
      RETURN
      END

5.5 代码分析流程图5.5 Code Analysis Flowchart

flowchart TD
    A[Input F_n, F_{n+1}, dt
输入变形梯度和时间步] --> B[Compute F_n^-1
计算旧变形梯度逆]
    B --> C[Compute L = (F_{n+1}-F_n)/dt * F_n^-1
计算速度梯度]
    C --> D[Extract W = 0.5*(L-L^T)
提取旋率张量]
    D --> E[Build A = I-0.5*W*dt
构造矩阵A]
    E --> F[Build B = I+0.5*W*dt
构造矩阵B]
    F --> G[Compute A^-1
计算A的逆]
    G --> H[DROT = A^-1 * B
计算旋转增量]
    H --> I[Output orthogonal DROT
输出正交旋转张量]

Chapter 6 核心本构算法Chapter 6 Core Constitutive Algorithm

核心本构算法CRYSTALPLASTICITY是黄永刚晶体塑性模型的计算核心。它接收单个积分点的应力、应变增量、状态变量和旋转张量，执行晶体塑性本构积分，返回更新后的应力和状态变量。The core constitutive algorithm CRYSTALPLASTICITY is the computational core of the Huang crystal plasticity model. It receives the stress, strain increment, state variables, and rotation tensor for a single integration point, performs crystal plasticity constitutive integration, and returns the updated stress and state variables.

6.1 CRYSTALPLASTICITY流程6.1 CRYSTALPLASTICITY Flow

flowchart TD
    A[Input: stress, dstrain, statev, DROT
输入应力应变状态旋转] --> B[Rotate to crystal frame
旋转到晶体坐标系]
    B --> C[Read Euler angles from STATEV
读取欧拉角]
    C --> D[Compute slip directions/normals
计算滑移方向/法向]
    D --> E[Compute Schmid tensors
计算Schmid张量]
    E --> F[Compute elastic stiffness [C]
计算弹性刚度矩阵]
    F --> G[Initialize trial stress
初始化试探应力]
    G --> H{Newton-Raphson iteration
牛顿迭代}
    H --> I[Compute resolved shear stress tau^alpha
计算分切应力]
    I --> J[Compute shear rate gammadot^alpha
计算剪切应变率]
    J --> K[Compute plastic strain rate
计算塑性应变率]
    K --> L[Compute stress residual
计算应力残差]
    L --> M{Converged? ||residual||<tol
收敛判断}
    M -->|No| N[Compute Jacobian and update
计算Jacobian并更新]
    N --> H
    M -->|Yes| O[Update state variables
更新状态变量]
    O --> P[Rotate back to global frame
旋转回全局坐标系]
    P --> Q[Output new stress & statev
输出新应力和状态]

6.2 弹性刚度矩阵6.2 Elastic Stiffness

晶体塑性模型中的弹性响应通常是各向异性的，由弹性刚度矩阵C描述。对于立方晶体（如FCC和BCC），弹性刚度矩阵在晶体坐标系中只有三个独立常数：C₁₁, C₁₂, C₄₄。The elastic response in crystal plasticity models is typically anisotropic, described by the elastic stiffness matrix C. For cubic crystals (such as FCC and BCC), the elastic stiffness matrix in the crystal coordinate system has only three independent constants: C₁₁, C₁₂, C₄₄.

      SUBROUTINE ElasticStiffness(C11,C12,C44,CC)
      DOUBLE PRECISION C11,C12,C44,CC(6,6)
      INTEGER I,J
      
      DO I=1,6
        DO J=1,6
          CC(I,J)=0.0
        END DO
      END DO
      
      ! Cubic crystal elastic stiffness in crystal coordinates
      ! 立方晶体弹性刚度（晶体坐标系）
      CC(1,1)=C11;  CC(2,2)=C11;  CC(3,3)=C11
      CC(1,2)=C12;  CC(1,3)=C12;  CC(2,3)=C12
      CC(2,1)=C12;  CC(3,1)=C12;  CC(3,2)=C12
      CC(4,4)=C44;  CC(5,5)=C44;  CC(6,6)=C44
      
      RETURN
      END

对于六方晶体（HCP），弹性刚度矩阵有五个独立常数：C₁₁, C₁₂, C₁₃, C₃₃, C₄₄。弹性刚度矩阵必须在晶体坐标系中定义，然后通过取向变换矩阵旋转到全局坐标系。For hexagonal crystals (HCP), the elastic stiffness matrix has five independent constants: C₁₁, C₁₂, C₁₃, C₃₃, C₄₄. The elastic stiffness matrix must be defined in the crystal coordinate system and then rotated to the global coordinate system through the orientation transformation matrix.

6.3 滑移系生成6.3 Slip System Generation

滑移系由滑移方向s和滑移面法向n定义。对于FCC晶体，有12个滑移系，属于{111}<110>滑移系族。对于BCC晶体，滑移系通常包括{110}<111>、{112}<111>和{123}<111>族。Slip systems are defined by slip direction s and slip plane normal n. For FCC crystals, there are 12 slip systems belonging to the {111}<110> slip system family. For BCC crystals, slip systems typically include {110}<111>, {112}<111>, and {123}<111> families.

      ! FCC crystal slip systems: 12 systems
      ! FCC晶体滑移系：12个
      ! Slip planes: {111} family
      ! Slip directions: <110> family
      
      ! Normalized slip plane normals (n)
      DATA SNDATA /
     &  0.70710678, 0.70710678, 0.0,
     &  0.70710678,-0.70710678, 0.0,
     &  0.70710678, 0.0, 0.70710678,
     &  0.70710678, 0.0,-0.70710678,
     &  0.0, 0.70710678, 0.70710678,
     &  0.0, 0.70710678,-0.70710678,
     & -0.70710678, 0.70710678, 0.0,
     & -0.70710678,-0.70710678, 0.0,
     & -0.70710678, 0.0, 0.70710678,
     & -0.70710678, 0.0,-0.70710678,
     &  0.0,-0.70710678, 0.70710678,
     &  0.0,-0.70710678,-0.70710678 /
      
      ! Normalized slip directions (s)
      DATA SLDATA /
     &  0.70710678,-0.70710678, 0.0,
     &  0.70710678, 0.70710678, 0.0,
     &  0.70710678, 0.0,-0.70710678,
     &  0.70710678, 0.0, 0.70710678,
     &  0.0, 0.70710678,-0.70710678,
     &  0.0, 0.70710678, 0.70710678,
     &  0.70710678, 0.70710678, 0.0,
     &  0.70710678,-0.70710678, 0.0,
     &  0.70710678, 0.0, 0.70710678,
     &  0.70710678, 0.0,-0.70710678,
     &  0.0, 0.70710678, 0.70710678,
     &  0.0, 0.70710678,-0.70710678 /

Schmid张量（取向因子矩阵）定义为P^α = 0.5(s^α⊗n^α + n^α⊗s^α)。分切应力τ^α = σ : P^α，即应力张量在Schmid张量上的投影。The Schmid tensor (orientation factor matrix) is defined as P^α = 0.5(s^α⊗n^α + n^α⊗s^α). The resolved shear stress τ^α = σ : P^α, which is the projection of the stress tensor onto the Schmid tensor.

6.4 隐式积分6.4 Implicit Integration

尽管VUMAT运行在显式求解器中，但每个积分点内部的本构积分通常采用隐式或半隐式格式，以保证数值稳定性。黄永刚模型使用广义中点法则（Theta方法）进行应力更新。Although VUMAT runs in an explicit solver, the constitutive integration within each integration point typically uses implicit or semi-implicit format to ensure numerical stability. The Huang model uses the generalized midpoint rule (Theta method) for stress update.

其中θ为积分参数，θ=1对应完全隐式（向后Euler），θ=0对应显式（向前Euler），θ=0.5对应中点法则。黄永刚原始代码通常使用θ=0.5，兼顾稳定性和精度。Where θ is the integration parameter: θ=1 corresponds to fully implicit (backward Euler), θ=0 corresponds to explicit (forward Euler), and θ=0.5 corresponds to the midpoint rule. The original Huang code typically uses θ=0.5, balancing stability and accuracy.

6.5 Newton-Raphson迭代6.5 Newton-Raphson Iteration

在每个增量步中，需要通过Newton-Raphson迭代求解非线性方程组，未知数是所有滑移系的剪切应变率γ̇^α。残差向量定义为应力更新方程的误差。In each increment, a system of nonlinear equations must be solved via Newton-Raphson iteration, where the unknowns are the shear strain rates γ̇^α on all slip systems. The residual vector is defined as the error in the stress update equation.

      ! Newton-Raphson iteration for stress update
      ! 应力更新的Newton-Raphson迭代
      
      ITER=0
      MAXITER=50
      TOLER=1.0E-8
      
      DO WHILE (ITER .LT. MAXITER)
        ITER=ITER+1
        
        ! Compute resolved shear stress on all slip systems
        ! 计算所有滑移系上的分切应力
        DO ISLIP=1,NSLIP
          TAU(ISLIP)=0.0
          DO I=1,6
            TAU(ISLIP)=TAU(ISLIP)+STRESS(I)*SCHMID(I,ISLIP)
          END DO
        END DO
        
        ! Compute shear strain rate from power law
        ! 从幂律计算剪切应变率
        DO ISLIP=1,NSLIP
          RATIO=ABS(TAU(ISLIP))/G(ISLIP)
          IF (RATIO .GT. 1.0E-10) THEN
            GAMMADOT(ISLIP)=GAM0*SIGN(RATIO**RATEEXP,TAU(ISLIP))
          ELSE
            GAMMADOT(ISLIP)=0.0
          END IF
        END DO
        
        ! Compute stress residual
        ! 计算应力残差
        DO I=1,6
          DEPLAS(I)=0.0
          DO ISLIP=1,NSLIP
            DEPLAS(I)=DEPLAS(I)+DT*GAMMADOT(ISLIP)*SCHMID(I,ISLIP)
          END DO
          STRESS_TRIAL(I)=STRESS_OLD(I)+DDSTRESS_ELAS(I)
          RESID(I)=STRESS(I)-STRESS_TRIAL(I)+DOT_CC_DEPLAS(I)
        END DO
        
        ! Check convergence
        ! 检查收敛
        RESNORM=0.0
        DO I=1,6
          RESNORM=RESNORM+RESID(I)**2
        END DO
        RESNORM=SQRT(RESNORM)
        
        IF (RESNORM .LT. TOLER) THEN
          CONVERGED=.TRUE.
          EXIT
        END IF
        
        ! Compute Jacobian and solve for correction
        ! 计算Jacobian并求解修正量
        CALL ComputeJacobian(...)
        CALL SolveLinearSystem(...)
        
        ! Update stress
        ! 更新应力
        DO I=1,6
          STRESS(I)=STRESS(I)+DELTA_STRESS(I)
        END DO
        
      END DO

Chapter 7 显式积分应力更新Chapter 7 Stress Update in Explicit Integration

在黄永刚VUMAT实现中，应力更新是核心计算环节。它结合了显式时间推进和隐式本构积分，在每个时间增量步内通过Newton-Raphson迭代确保本构方程的满足。In the Huang VUMAT implementation, stress update is the core computational step. It combines explicit time advancement with implicit constitutive integration, ensuring constitutive equation satisfaction within each time increment through Newton-Raphson iteration.

7.1 幂律本构7.1 Power Law Constitutive

滑移系上的剪切应变率通常采用幂律（Power Law）或双曲正弦（Hyperbolic Sine）形式描述。幂律形式最为常用，其表达式为：The shear strain rate on slip systems is typically described by a power law or hyperbolic sine form. The power law form is most commonly used, expressed as:

其中γ̇₀是参考剪切应变率，g^α是滑移阻力（临界分切应力），m是应变率敏感性指数。当m趋近于0时，材料表现为率无关；m越大，率敏感性越强。Where γ̇_0 is the reference shear strain rate, g^α is the slip resistance (critical resolved shear stress), and m is the strain rate sensitivity exponent. As m approaches 0, the material behaves in a rate-independent manner; the larger m is, the stronger the rate sensitivity.

      ! Power law constitutive relation for slip rate
      ! 滑移应变率的幂律本构关系
      
      DO ISLIP=1,NSLIP
        ! Resolved shear stress / 分切应力
        TAU = STRESS(1)*SCHMID(1,ISLIP)
     &      + STRESS(2)*SCHMID(2,ISLIP)
     &      + STRESS(3)*SCHMID(3,ISLIP)
     &      + STRESS(4)*SCHMID(4,ISLIP)
     &      + STRESS(5)*SCHMID(5,ISLIP)
     &      + STRESS(6)*SCHMID(6,ISLIP)
        
        ! Slip resistance from state variable / 滑移阻力
        G_SLIP = STATEV(45+ISLIP)
        
        ! Power law / 幂律
        IF (ABS(G_SLIP) .GT. 1.0E-12) THEN
          RATIO = ABS(TAU) / G_SLIP
          GAMMADOT(ISLIP) = GAM0 * (RATIO**EXPON)
     &                    * SIGN(1.0, TAU)
        ELSE
          GAMMADOT(ISLIP) = 0.0
        END IF
      END DO

7.2 Theta隐式积分7.2 Theta Implicit Scheme

应力更新采用广义中点法则。将应变率分解为弹性和塑性部分，弹性部分通过Hooke定律与应力率关联，塑性部分由所有激活滑移系的剪切贡献叠加。Stress update uses the generalized midpoint rule. The strain rate is decomposed into elastic and plastic parts; the elastic part is related to the stress rate through Hooke’s law, and the plastic part is the superposition of shear contributions from all active slip systems.

当θ=0.5时，该格式具有二阶时间精度，且无条件稳定（在本构积分层面）。由于显式求解器整体是条件稳定的，本构积分的稳定性不影响整体稳定性，但会影响每个积分点内部计算的准确性。When θ=0.5, this scheme has second-order temporal accuracy and is unconditionally stable (at the constitutive integration level). Since the explicit solver as a whole is conditionally stable, the stability of constitutive integration does not affect overall stability but does affect the accuracy of internal calculations at each integration point.

7.3 线性化求解7.3 Linearized Solution

Newton-Raphson迭代的核心是构建并求解线性化方程组。Jacobian矩阵由应力对剪切应变率的导数构成。对于幂律本构，Jacobian可以解析求导，避免了数值差分的精度损失和额外计算开销。The core of Newton-Raphson iteration is constructing and solving the linearized equation system. The Jacobian matrix consists of derivatives of stress with respect to shear strain rates. For the power law constitutive relation, the Jacobian can be derived analytically, avoiding accuracy loss and additional computational overhead from numerical differentiation.

      ! Analytical Jacobian for power-law crystal plasticity
      ! 幂律晶体塑性的解析Jacobian
      !
      ! d(RESID)/d(GAMMADOT) = dt * C : P^alpha : P^beta
      !     + dt * (dt * theta) * (dgammadot/dtau) * (C:P^alpha) * P^beta
      
      DO ISLIP=1,NSLIP
        DO JSLIP=1,NSLIP
          JACOB(ISLIP,JSLIP)=0.0
          
          ! Elastic contribution
          DO I=1,6
            DO J=1,6
              JACOB(ISLIP,JSLIP)=JACOB(ISLIP,JSLIP)
     &          +DT*CC(I,J)*SCHMID(J,JSLIP)*SCHMID(I,ISLIP)
            END DO
          END DO
          
          ! Plastic contribution (rate sensitivity)
          TAU_J=0.0
          DO I=1,6
            TAU_J=TAU_J+STRESS(I)*SCHMID(I,JSLIP)
          END DO
          
          IF (ABS(TAU_J) .GT. 1.0E-12) THEN
            DGDTAU=GAM0/EXPON
     &            *(ABS(TAU_J)/G(JSLIP))**(1.0/EXPON-1.0)
     &            /G(JSLIP)
            
            DO I=1,6
              CP=0.0
              DO J=1,6
                CP=CP+CC(I,J)*SCHMID(J,JSLIP)
              END DO
              JACOB(ISLIP,JSLIP)=JACOB(ISLIP,JSLIP)
     &          +DT*DT*THETA*DGDTAU*CP*SCHMID(I,ISLIP)
     &          *SIGN(1.0,TAU_J)
            END DO
          END IF
          
        END DO
      END DO

7.4 应力更新公式7.4 Stress Update Formula

完整的应力更新流程可以总结为以下步骤：首先计算弹性试探应力，然后通过迭代调整应力使其满足本构关系，最后将更新后的应力存储到STATEV并旋转回全局坐标系。The complete stress update procedure can be summarized in the following steps: first compute the elastic trial stress, then iteratively adjust the stress to satisfy the constitutive relation, and finally store the updated stress in STATEV and rotate back to the global coordinate system.

上述公式体现了晶体塑性应力更新的本质：从弹性试探应力中减去塑性应变对应的应力松弛量。由于多个滑移系可能同时激活，各滑移系的剪切应变率通过Newton-Raphson迭代耦合求解。The above formula reflects the essence of crystal plasticity stress update: subtracting the stress relaxation corresponding to plastic strain from the elastic trial stress. Since multiple slip systems may be simultaneously active, the shear strain rates of each slip system are solved in a coupled manner through Newton-Raphson iteration.

Chapter 8 材料参数Chapter 8 Material Parameters

黄永刚VUMAT使用的材料参数与UMAT版本完全相同。这些参数通过PROPS数组传入子程序，在Abaqus输入文件中使用*USER MATERIAL关键字定义。Huang’s VUMAT uses exactly the same material parameters as the UMAT version. These parameters are passed into the subroutine through the PROPS array and are defined in the Abaqus input file using the *USER MATERIAL keyword.

参数表中前几个参数通常用于定义晶体类型和弹性常数，后续参数定义塑性流动和硬化行为。理解每个参数的物理意义对于正确模拟材料行为至关重要。The first few parameters in the table typically define the crystal type and elastic constants, while subsequent parameters define plastic flow and hardening behavior. Understanding the physical meaning of each parameter is crucial for correctly simulating material behavior.

硬化演化采用Voce硬化律描述，滑移阻力g^α随累积剪切应变γ的演化关系为：Hardening evolution is described by the Voce hardening law. The evolution of slip resistance g^α with accumulated shear strain γ is:

其中γ = Σ_α |γ^α| 是所有滑移系上累积剪切应变的总和。自硬化和潜硬化通过硬化矩阵h^αβ描述，对角元素为自硬化，非对角元素为潜硬化，通常取q·h_self。Where γ = Σ_α |γ^α| is the total accumulated shear strain over all slip systems. Self-hardening and latent hardening are described by the hardening matrix h^αβ, where diagonal elements are self-hardening and off-diagonal elements are latent hardening, typically taken as q·h_self.

Chapter 9 关键代码Chapter 9 Key Code

本章展示黄永刚VUMAT的关键代码片段，包括外壳代码、DROT计算和应变转换。代码注释采用中英双语，以便不同语言背景的读者理解。This chapter presents key code snippets of the Huang VUMAT, including the wrapper code, DROT computation, and strain conversion. Code comments are bilingual to facilitate readers from different language backgrounds.

9.1 VUMAT外壳代码9.1 VUMAT Wrapper Code

      SUBROUTINE VUMAT(
     &  NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, LANNEAL,
     &  STEPTIME, TOTALTIME, DT, CMNAME, COORDMP, CHARLENGTH,
     &  PROPS, DENSITY, STRAININC, RELSPININC,
     &  TEMPOLD, STRETCHOLD, DEFGRADOLD, FIELDOLD,
     &  STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,
     &  TEMPNEW,  STRETCHNEW,  DEFGRADNEW,  FIELDNEW,
     &  STRESSNEW,  STATENEW,  ENERINTERNNEW, ENERINELASNEW )
      
      ! ================================================================
      ! VUMAT wrapper for crystal plasticity (Huang's model)
      ! 黄永刚晶体塑性VUMAT外壳子程序
      ! ================================================================
      
      IMPLICIT NONE
      
      ! Arguments / 参数声明
      INTEGER NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, LANNEAL
      DOUBLE PRECISION STEPTIME, TOTALTIME, DT
      CHARACTER*80 CMNAME
      DOUBLE PRECISION PROPS(NPROPS), DENSITY(NBLOCK)
      DOUBLE PRECISION COORDMP(NBLOCK), CHARLENGTH(NBLOCK)
      DOUBLE PRECISION STRAININC(NBLOCK,NDIR+NSHR)
      DOUBLE PRECISION RELSPININC(NBLOCK,NSHR)
      DOUBLE PRECISION TEMPOLD(NBLOCK), STRETCHOLD(NBLOCK,NDIR+NSHR)
      DOUBLE PRECISION DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR)
      DOUBLE PRECISION FIELDOLD(NBLOCK,NFIELDV)
      DOUBLE PRECISION STRESSOLD(NBLOCK,NDIR+NSHR)
      DOUBLE PRECISION STATEOLD(NBLOCK,NSTATEV)
      DOUBLE PRECISION ENERINTERNOLD(NBLOCK), ENERINELASOLD(NBLOCK)
      DOUBLE PRECISION TEMPNEW(NBLOCK), STRETCHNEW(NBLOCK,NDIR+NSHR)
      DOUBLE PRECISION DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR)
      DOUBLE PRECISION FIELDNEW(NBLOCK,NFIELDV)
      DOUBLE PRECISION STRESSNEW(NBLOCK,NDIR+NSHR)
      DOUBLE PRECISION STATENEW(NBLOCK,NSTATEV)
      DOUBLE PRECISION ENERINTERNNEW(NBLOCK), ENERINELASNEW(NBLOCK)
      
      ! Local variables / 局部变量
      INTEGER KM, K1, K2, K3
      DOUBLE PRECISION STRESS(6), DSTRAN(6), STATEV(1000)
      DOUBLE PRECISION FOLD(3,3), FNEW(3,3), DROT(3,3)
      DOUBLE PRECISION PROPS_LOC(50)
      DOUBLE PRECISION TEMPK
      
      ! Copy material properties / 复制材料参数
      DO K1=1,NPROPS
        PROPS_LOC(K1)=PROPS(K1)
      END DO
      
      ! ================================================================
      ! Loop over all integration points in the block
      ! 遍历块中所有积分点
      ! ================================================================
      DO KM=1,NBLOCK
        
        ! Extract old stress / 提取旧应力
        DO K1=1,NDIR+NSHR
          STRESS(K1)=STRESSOLD(KM,K1)
        END DO
        
        ! Extract old state variables / 提取旧状态变量
        DO K1=1,NSTATEV
          STATEV(K1)=STATEOLD(KM,K1)
        END DO
        
        ! Extract deformation gradients / 提取变形梯度
        CALL EXTRACT_F(DEFGRADOLD,KM,FOLD,NDIR,NSHR)
        CALL EXTRACT_F(DEFGRADNEW,KM,FNEW,NDIR,NSHR)
        
        ! Compute rotation tensor DROT / 计算旋转张量
        CALL COMPUTEDROT(FOLD,FNEW,DT,DROT)
        
        ! Convert strain increment / 转换应变增量
        DO K1=1,NDIR+NSHR
          DSTRAN(K1)=STRAININC(KM,K1)
        END DO
        ! Tensor shear strain conversion / 张量剪切应变转换
        DO K1=NDIR+1,NDIR+NSHR
          DSTRAN(K1)=DSTRAN(K1)*0.5
        END DO
        
        ! Temperature / 温度
        TEMPK=TEMPOLD(KM)
        
        ! Call core constitutive subroutine / 调用核心本构子程序
        CALL CRYSTALPLASTICITY(
     &    STRESS,DSTRAN,STATEV,NSTATEV,PROPS_LOC,NPROPS,
     &    DROT,NDIR,NSHR,DT,TEMPK)
        
        ! Store updated stress / 存储更新后的应力
        DO K1=1,NDIR+NSHR
          STRESSNEW(KM,K1)=STRESS(K1)
        END DO
        
        ! Store updated state variables / 存储更新后的状态变量
        DO K1=1,NSTATEV
          STATENEW(KM,K1)=STATEV(K1)
        END DO
        
      END DO
      
      RETURN
      END

9.2 DROT计算代码9.2 DROT Code

      SUBROUTINE COMPUTEDROT(FN,FN1,DT,DROT)
      ! ================================================================
      ! Compute incremental rotation tensor from deformation gradient
      ! 从变形梯度计算增量旋转张量
      ! Using 2nd-order Cayley transform
      ! 使用二阶Cayley变换
      ! ================================================================
      IMPLICIT NONE
      DOUBLE PRECISION FN(3,3),FN1(3,3),DT,DROT(3,3)
      DOUBLE PRECISION FINV(3,3),L(3,3),W(3,3)
      DOUBLE PRECISION A(3,3),B(3,3),DET,WVEC(3)
      INTEGER I,J,K
      
      ! Compute inverse of FN / 计算FN的逆
      CALL M3INV(FN,FINV,DET)
      IF (ABS(DET) .LT. 1.0E-12) THEN
        WRITE(*,*) 'ERROR: Singular deformation gradient'
        CALL XIT
      END IF
      
      ! Velocity gradient: L = dF/dt * F^-1 / 速度梯度
      DO I=1,3
        DO J=1,3
          L(I,J)=0.0
          DO K=1,3
            L(I,J)=L(I,J)+(FN1(I,K)-FN(I,K))/DT*FINV(K,J)
          END DO
        END DO
      END DO
      
      ! Spin tensor W = 0.5*(L - L^T) / 旋率张量
      W(1,1)=0.0
      W(2,2)=0.0
      W(3,3)=0.0
      W(1,2)=0.5*(L(1,2)-L(2,1))
      W(1,3)=0.5*(L(1,3)-L(3,1))
      W(2,3)=0.5*(L(2,3)-L(3,2))
      W(2,1)=-W(1,2)
      W(3,1)=-W(1,3)
      W(3,2)=-W(2,3)
      
      ! Extract axial vector / 提取轴矢量
      WVEC(1)=W(3,2)
      WVEC(2)=W(1,3)
      WVEC(3)=W(2,1)
      
      ! Build Cayley transform matrices / 构造Cayley变换矩阵
      DO I=1,3
        DO J=1,3
          A(I,J)=-0.5*DT*W(I,J)
          B(I,J)= 0.5*DT*W(I,J)
        END DO
        A(I,I)=A(I,I)+1.0
        B(I,I)=B(I,I)+1.0
      END DO
      
      ! Invert A / 求A的逆
      CALL M3INV(A,FINV,DET)
      
      ! DROT = A^-1 * B / 计算旋转增量
      DO I=1,3
        DO J=1,3
          DROT(I,J)=0.0
          DO K=1,3
            DROT(I,J)=DROT(I,J)+FINV(I,K)*B(K,J)
          END DO
        END DO
      END DO
      
      RETURN
      END

9.3 应变转换代码9.3 Strain Conversion Code

      SUBROUTINE CONVERT_STRAIN(STRAININC,KM,NDIR,NSHR,DSTRAN)
      ! ================================================================
      ! Convert strain increment from Abaqus/Explicit format to
      ! tensor shear strain format required by constitutive algorithm
      ! 将Abaqus/Explicit的应变增量格式转换为核心本构算法需要的张量剪切应变格式
      ! ================================================================
      IMPLICIT NONE
      INTEGER KM, NDIR, NSHR
      DOUBLE PRECISION STRAININC(NBLOCK,NDIR+NSHR), DSTRAN(6)
      INTEGER K1
      
      ! Copy direct strains / 复制正应变分量
      DO K1=1,NDIR
        DSTRAN(K1)=STRAININC(KM,K1)
      END DO
      
      ! Convert engineering shear to tensor shear
      ! 工程剪切应变 -> 张量剪切应变
      ! Abaqus/Explicit: engineering shear strain = 2 * tensor shear
      ! 工程剪切应变 = 2 * 张量剪切应变
      DO K1=NDIR+1,NDIR+NSHR
        DSTRAN(K1)=STRAININC(KM,K1)*0.5
      END DO
      
      RETURN
      END
      
      ! ================================================================
      ! Reverse conversion for output (if needed)
      ! 输出时的反向转换（如需要）
      ! ================================================================
      SUBROUTINE CONVERT_STRAIN_BACK(DSTRAN,NDIR,NSHR,STRAINOUT)
      IMPLICIT NONE
      INTEGER NDIR, NSHR
      DOUBLE PRECISION DSTRAN(6), STRAINOUT(6)
      INTEGER K1
      
      DO K1=1,NDIR
        STRAINOUT(K1)=DSTRAN(K1)
      END DO
      DO K1=NDIR+1,NDIR+NSHR
        STRAINOUT(K1)=DSTRAN(K1)*2.0
      END DO
      
      RETURN
      END

Chapter 10 使用注意事项Chapter 10 Usage Notes

成功使用黄永刚VUMAT子程序需要注意以下事项。这些建议来源于实际工程应用经验和常见问题排查。Successfully using the Huang VUMAT subroutine requires attention to the following items. These recommendations come from practical engineering application experience and common troubleshooting.

10.1 输入文件设置10.1 Input File Setup

• 必须在*MATERIAL段中使用*USER MATERIAL, TYPE=MECHANICAL, CONSTANTS=N定义材料参数Use *USER MATERIAL, TYPE=MECHANICAL, CONSTANTS=N in the *MATERIAL section
• DEPVAR必须大于等于STATEV数组实际使用的长度DEPVAR must be greater than or equal to the actual length used by the STATEV array
• 使用*ORIENTATION定义晶体取向，或在STATEV中存储欧拉角Define crystal orientation using *ORIENTATION, or store Euler angles in STATEV
• 单元类型应选用显式兼容类型，如C3D8RUse explicit-compatible element types such as C3D8R
• 必须在*STEP中使用*DYNAMIC, EXPLICITUse *DYNAMIC, EXPLICIT in the *STEP section

10.2 常见错误与排查10.2 Common Errors and Troubleshooting

• 错误：Signal 11 (Segmentation fault) — 通常由STATEV越界访问引起。检查DEPVAR和数组维度。Error: Signal 11 — usually caused by out-of-bounds STATEV access. Check DEPVAR and array dimensions.
• 错误：Zero pivot in DROT计算 — 变形梯度可能奇异（体积压缩至零）。检查网格畸变。Error: Zero pivot in DROT — deformation gradient may be singular. Check mesh distortion.
• 错误：应力不收敛 — 时间步长过大或材料参数不合理。减小时间步或调整率敏感指数m。Error: Stress not converging — time increment too large or material parameters unreasonable. Reduce time increment or adjust rate sensitivity m.
• 错误：织构演化异常 — DROT计算可能有误。验证旋转张量的正交性。Error: Abnormal texture evolution — DROT computation may be incorrect. Verify orthogonality of rotation tensor.
• 错误：剪切应力为预期值的两倍 — 忘记进行剪切应变转换。检查DSTRAN(K1)=STRAININC(KM,K1)*0.5。Error: Shear stress twice expected — forgot shear strain conversion. Check DSTRAN=STRAININC*0.5.

10.3 性能优化建议10.3 Performance Optimization

• 使用编译器优化选项（-O3）编译子程序Compile subroutine with optimizer flags (-O3)
• 尽量减少子程序内部的动态内存分配Minimize dynamic memory allocation inside subroutine
• 对于大规模模拟，考虑使用并行版本（VUMAT支持Abaqus自动向量化）For large-scale simulations, consider parallel versions
• 状态变量数目影响内存占用，不必要的变量不要存入STATEVState variable count affects memory; avoid storing unnecessary variables
• 如果不需要所有滑移系的历史，可以减少输出变量Reduce output variables if not all slip system histories are needed

10.4 验证与确认10.4 Verification and Validation

• 单晶单轴拉伸：验证应力-应变曲线与实验或解析解一致Single-crystal uniaxial tension: verify stress-strain curve matches experiments or analytical solutions
• 简单剪切：验证剪切应力-剪切应变关系和取向转动Simple shear: verify shear stress-strain relation and orientation rotation
• 多晶泰勒模型：与均匀化理论结果对比Polycrystal Taylor model: compare with homogenization theory results
• 网格收敛性：不同网格密度下的结果应趋于一致Mesh convergence: results should converge with refined mesh
• 能量守恒：弹性卸载后应变能应恢复Energy conservation: elastic strain energy should recover upon unloading

10.5 学术引用10.5 Academic Citation

使用黄永刚晶体塑性模型时，建议在论文中引用以下原始文献：When using the Huang crystal plasticity model, it is recommended to cite the following original literature:

• Huang Y. A user-material subroutine incorporating single crystal plasticity in the Abaqus finite element program. Harvard University Report, 1991.Huang Y. A user-material subroutine incorporating single crystal plasticity in the Abaqus finite element program. Harvard University Report, 1991.
• Peirce D, Asaro RJ, Needleman A. An analysis of nonuniform and localized deformation in ductile single crystals. Acta Metallurgica, 1982, 30(6): 1087-1119.Peirce D, Asaro RJ, Needleman A. An analysis of nonuniform and localized deformation in ductile single crystals. Acta Metallurgica, 1982, 30(6): 1087-1119.
• Asaro RJ. Crystal plasticity. Journal of Applied Mechanics, 1983, 50(4b): 921-934.Asaro RJ. Crystal plasticity. Journal of Applied Mechanics, 1983, 50(4b): 921-934.

附录A 数学符号表Appendix A Mathematical Symbols

本节列出本文档中使用的所有数学符号，便于读者快速查阅。This section lists all mathematical symbols used in this document for quick reference.

附录B Fortran编译指南Appendix B Fortran Compilation Guide

B.1 编译器选择B.1 Compiler Selection

Abaqus/Explicit支持多种Fortran编译器。对于Abaqus 2020及以前版本，通常使用Intel Fortran Compiler (ifort)。对于Abaqus 2021及以后版本，推荐使用Intel oneAPI中的ifx编译器。Abaqus/Explicit supports multiple Fortran compilers. For Abaqus 2020 and earlier, the Intel Fortran Compiler (ifort) is typically used. For Abaqus 2021 and later, the ifx compiler from Intel oneAPI is recommended.

B.2 编译命令B.2 Compilation Commands

在Windows平台上，可以使用Abaqus自带的编译环境或手动编译。以下是常见的编译命令示例：On Windows, you can use the Abaqus-provided compilation environment or compile manually. The following are common compilation command examples:

! 在Abaqus Command窗口中执行
! Execute in Abaqus Command window
abaqus make library=vumat.f90

! 或者编译为对象文件然后链接
! Or compile to object file then link
ifort /c /O3 /Qip vumat.f90
link /DLL vumat.obj

! 64位编译
! 64-bit compilation
ifort /c /O3 /Qip /heap-arrays /QaxCORE-AVX2 vumat.f90

! 生成动态链接库
! Generate dynamic link library
link /DLL /OUT:vumat.dll vumat.obj

B.3 编译优化选项B.3 Compiler Optimization Flags

• /O3 — 最高级别优化 / Maximum optimization level/O3 — Maximum optimization level
• /Qip — 过程间优化 / Interprocedural optimization/Qip — Interprocedural optimization
• /heap-arrays — 大数组自动分配堆内存 / Auto heap allocation for large arrays/heap-arrays — Auto heap allocation for large arrays
• /QaxCORE-AVX2 — 启用AVX2向量化 / Enable AVX2 vectorization/QaxCORE-AVX2 — Enable AVX2 vectorization
• /fp:precise — 保证浮点精度 / Ensure floating-point precision/fp:precise — Ensure floating-point precision

附录C Abaqus输入文件示例Appendix C Abaqus Input File Example

以下是一个完整的Abaqus/Explicit输入文件片段，展示如何定义使用黄永刚VUMAT子程序的材料和单元。The following is a complete Abaqus/Explicit input file snippet showing how to define materials and elements using the Huang VUMAT subroutine.

*Heading
** Crystal plasticity explicit analysis example
** 晶体塑性显式分析示例
**
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** Part definition
*Part, name=Specimen
*Node
      1,  0.0,  0.0,  0.0
      2,  1.0,  0.0,  0.0
      3,  1.0,  1.0,  0.0
      4,  0.0,  1.0,  0.0
      5,  0.0,  0.0,  1.0
      6,  1.0,  0.0,  1.0
      7,  1.0,  1.0,  1.0
      8,  0.0,  1.0,  1.0
*Element, type=C3D8R
 1, 1, 2, 3, 4, 5, 6, 7, 8
*End Part
**
** Material definition with VUMAT
** 使用VUMAT的材料定义
*Material, name=Aluminum-Single-Crystal
*Density
 2.70e-9
*User Material, type=MECHANICAL, constants=11
** PROPS(1): Crystal type = 1 (FCC)
** PROPS(2-4): C11, C12, C44 (MPa)
** PROPS(5): gamma_dot_0 (1/s)
** PROPS(6): rate sensitivity m
** PROPS(7): initial slip resistance g0 (MPa)
** PROPS(8): saturation slip resistance gs (MPa)
** PROPS(9): hardening coefficient h0 (MPa)
** PROPS(10): hardening exponent a
** PROPS(11): latent hardening q
 1.0, 106750.0, 60410.0, 28340.0, 1.0, 0.02, 31.0, 63.0, 75.0, 1.3, 1.0
*Depvar
 70
**
** Section assignment
*Solid Section, elset=AllElements, material=Aluminum-Single-Crystal
 1.0
*End Instance
**
** Boundary conditions
*Boundary
 Set-Bottom, ENCASTRE
**
** Step definition
*Step, name=Loading, nlgeom=YES
*Dynamic, Explicit
, 1.0e-6
** 显式动力学分析步，时间步长自动计算
** Explicit dynamic step, automatic time increment
*Bulk Viscosity
 0.06, 1.2
**
** Boundary condition application
*Boundary, type=VELOCITY
 Set-Top, 2, 2, 1000.0
**
** Output requests
*Output, field, variable=PRESELECT
*Output, history, variable=PRESELECT
*Element Output, elset=AllElements
 SDV, S, E, PE
*End Step

附录D 常见问题FAQAppendix D Frequently Asked Questions

D.1 VUMAT和UMAT可以同时编译使用吗？D.1 Can VUMAT and UMAT be compiled and used simultaneously?

可以。只要在输入文件中分别为Standard和Explicit分析步指定不同的材料名，就可以在同一个模型中同时使用UMAT和VUMAT。但注意两者不能用于同一个分析步。Yes. As long as different material names are specified for Standard and Explicit steps in the input file, UMAT and VUMAT can be used simultaneously in the same model. However, they cannot be used in the same analysis step.

D.2 如何提取特定滑移系的剪切应变？D.2 How to extract shear strain on a specific slip system?

在Abaqus/Viewer中，使用Results -> Field Output，选择SDV变量。根据STATEV布局定义（参见3.4节），第10-21个SDV分量对应12个滑移系的累积剪切应变。In Abaqus/Viewer, use Results -> Field Output and select SDV variables. According to the STATEV layout definition (see Section 3.4), SDV components 10-21 correspond to accumulated shear strains on the 12 slip systems.

D.3 显式分析的时间步长如何控制？D.3 How to control time increment in explicit analysis?

Abaqus/Explicit会自动计算稳定时间增量。用户可以通过*FIXED TIME INCREMENTATION或调整质量缩放（mass scaling）来间接控制。不建议直接固定非常小的时间步长，因为这会导致计算成本过高。Abaqus/Explicit automatically calculates the stable time increment. Users can indirectly control it through *FIXED TIME INCREMENTATION or by adjusting mass scaling. Directly fixing a very small time increment is not recommended as it leads to excessive computational cost.

D.4 晶体取向如何随机分布？D.4 How to distribute crystal orientations randomly?

可以使用Python脚本生成随机欧拉角，然后通过*Orientation关键字或直接在STATEV初始条件中指定。对于多晶模拟，通常使用Dream.3D、Neper等软件生成代表性体积元（RVE），然后映射到有限元网格。Random Euler angles can be generated using Python scripts and specified through the *Orientation keyword or directly in STATEV initial conditions. For polycrystal simulations, software such as Dream.3D or Neper is typically used to generate representative volume elements (RVE), which are then mapped to the finite element mesh.

D.5 模拟结果与实验差异大怎么办？D.5 What if simulation results differ significantly from experiments?

首先检查材料参数是否合理，特别是率敏感指数m和硬化参数。其次检查网格密度是否足够（建议至少3-4个单元穿过晶粒）。最后检查边界条件是否正确施加。建议从单晶简单加载开始验证，逐步过渡到多晶复杂加载。First check whether material parameters are reasonable, especially the rate sensitivity exponent m and hardening parameters. Second, check whether mesh density is sufficient (recommend at least 3-4 elements across a grain). Finally, check whether boundary conditions are correctly applied. It is recommended to start verification with single-crystal simple loading and gradually transition to polycrystal complex loading.

D.6 VUMAT支持热-力耦合分析吗？D.6 Does VUMAT support thermo-mechanical coupled analysis?

VUMAT接口提供了TEMPOLD和TEMPNEW温度变量，因此可以在子程序中实现热-力耦合本构。但需要注意，显式求解器中的绝热加热效应需要通过状态变量或用户子程序自行计算和传递。The VUMAT interface provides TEMPOLD and TEMPNEW temperature variables, so thermo-mechanical coupled constitutive relations can be implemented in the subroutine. However, note that adiabatic heating effects in the explicit solver need to be calculated and passed through state variables or user subroutines.

D.7 子程序运行速度很慢怎么办？D.7 What if the subroutine runs very slowly?

检查以下几点：（1）Newton-Raphson迭代次数是否过多，可考虑增大收敛容差；（2）状态变量数目是否过多，减少不必要的输出；（3）是否启用了编译器优化；（4）是否使用了合适的单元类型（C3D8R通常比C3D8快）；（5）考虑使用并行计算（MPI或线程并行）。Check the following: (1) whether Newton-Raphson iterations are excessive, consider increasing convergence tolerance; (2) whether state variable count is too large, reduce unnecessary output; (3) whether compiler optimization is enabled; (4) whether appropriate element types are used (C3D8R is typically faster than C3D8); (5) consider using parallel computing (MPI or thread parallelism).

附录E 扩展资源与参考文献Appendix E Extended Resources and References

E.1 核心参考文献E.1 Core References

• Asaro RJ, Rice JR. Strain localization in ductile single crystals. Journal of the Mechanics and Physics of Solids, 1977, 25(5): 309-338.Asaro RJ, Rice JR. Strain localization in ductile single crystals. Journal of the Mechanics and Physics of Solids, 1977, 25(5): 309-338.
• Hutchinson JW. Bounds and self-consistent estimates for creep of polycrystalline materials. Proceedings of the Royal Society A, 1976, 348(1652): 101-127.Hutchinson JW. Bounds and self-consistent estimates for creep of polycrystalline materials. Proceedings of the Royal Society A, 1976, 348(1652): 101-127.
• Kalidindi SR, Bronkhorst CA, Anand L. Crystallographic texture evolution in bulk deformation processing of FCC metals. Journal of the Mechanics and Physics of Solids, 1992, 40(3): 537-569.Kalidindi SR, Bronkhorst CA, Anand L. Crystallographic texture evolution in bulk deformation processing of FCC metals. Journal of the Mechanics and Physics of Solids, 1992, 40(3): 537-569.
• Roters F, Eisenlohr P, Hantcherli L, et al. Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications. Acta Materialia, 2010, 58(4): 1152-1211.Roters F, Eisenlohr P, Hantcherli L, et al. Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications. Acta Materialia, 2010, 58(4): 1152-1211.
• Van Houtte P, Li S, Seefeldt M, Delannay L. Deformation texture prediction: from the Taylor model to the advanced self-consistent model. Acta Materialia, 2005, 53(11): 3225-3236.Van Houtte P, Li S, Seefeldt M, Delannay L. Deformation texture prediction: from the Taylor model to the advanced self-consistent model. Acta Materialia, 2005, 53(11): 3225-3236.

E.2 推荐学习资源E.2 Recommended Learning Resources

• DAMASK（Dusseldorf Advanced Material Simulation Kit）— 开源晶体塑性有限元软件，提供丰富的本构模型和前后处理工具。DAMASK (Dusseldorf Advanced Material Simulation Kit) — Open-source crystal plasticity finite element software providing rich constitutive models and pre/post-processing tools.
• MTEX — MATLAB工具箱，用于晶体学织构分析和取向数据处理。MTEX — MATLAB toolbox for crystallographic texture analysis and orientation data processing.
• Dream.3D — 开源微结构重建和网格生成软件，常用于生成多晶RVE。Dream.3D — Open-source microstructure reconstruction and mesh generation software, commonly used for generating polycrystal RVEs.
• Neper — 多晶微结构生成和网格划分工具，支持Voronoi tessellation。Neper — Polycrystal microstructure generation and meshing tool supporting Voronoi tessellation.
• Abaqus Documentation — Abaqus用户子程序参考手册（User Subroutines Reference Manual）中的VUMAT章节。Abaqus Documentation — VUMAT chapter in Abaqus User Subroutines Reference Manual.

E.3 相关会议与期刊E.3 Relevant Conferences and Journals

• 期刊：International Journal of Plasticity, Journal of the Mechanics and Physics of Solids, Acta Materialia, Computational MechanicsJournals: International Journal of Plasticity, Journal of the Mechanics and Physics of Solids, Acta Materialia, Computational Mechanics
• 会议：ICHSIP（International Conference on the Heterogeneity of Plastic Deformation）, Plasticity, TMS Annual MeetingConferences: ICHSIP (International Conference on the Heterogeneity of Plastic Deformation), Plasticity, TMS Annual Meeting

E.4 版本历史E.4 Version History

附录F 术语表Appendix F Glossary

补充技术说明Supplementary Technical Notes

S.1 显式与隐式积分对比S.1 Explicit vs Implicit Integration Comparison

在晶体塑性有限元方法中，时间积分格式的选择直接影响计算效率和稳定性。显式积分不需要形成和求解全局刚度矩阵，每个增量步的计算成本较低，但受CFL稳定性条件限制，时间步长必须小于临界值。隐式积分允许更大的时间步长，但需要迭代求解非线性方程组，且需要一致切线刚度矩阵。In crystal plasticity finite element methods, the choice of time integration format directly affects computational efficiency and stability. Explicit integration does not require forming and solving the global stiffness matrix, resulting in lower computational cost per increment, but is subject to the CFL stability condition where the time increment must be less than a critical value. Implicit integration allows larger time increments but requires iterative solution of nonlinear equation systems and a consistent tangent stiffness matrix.

VUMAT采用了一种混合策略：在全局层面使用显式时间推进，避免全局迭代；在局部（积分点）层面使用隐式或半隐式本构积分，保证应力更新的准确性和稳定性。这种策略兼顾了显式求解器的计算效率和隐式本构积分的数值稳定性，是晶体塑性显式分析的最佳实践。VUMAT adopts a hybrid strategy: explicit time advancement at the global level avoids global iteration; implicit or semi-implicit constitutive integration at the local (integration point) level ensures accuracy and stability of stress update. This strategy balances the computational efficiency of explicit solvers with the numerical stability of implicit constitutive integration, representing best practice for crystal plasticity explicit analysis.

S.2 客观性要求与Jaumann率S.2 Objectivity Requirement and Jaumann Rate

在大变形分析中，应力更新必须满足客观性（objectivity）或框架无差异性（frame indifference）要求。这意味着本构关系在任意刚体旋转下应保持不变。在晶体塑性中，客观性通过将应力率定义在刚体旋转的共旋坐标系中来实现。In large deformation analysis, stress update must satisfy objectivity or frame indifference requirements. This means constitutive relations should remain invariant under arbitrary rigid body rotations. In crystal plasticity, objectivity is achieved by defining stress rates in a co-rotational coordinate system attached to the rigid body rotation.

Jaumann应力率是最常用的客观应力率之一，定义为应力张量对时间的物质导数减去旋率W与应力的对易子。在黄永刚VUMAT中，客观性通过增量旋转张量DROT实现：每次本构积分前，将应力从全局坐标系旋转到无旋坐标系；积分完成后，再旋转回全局坐标系。这种方法等价于使用Jaumann率进行应力更新。The Jaumann stress rate is one of the most commonly used objective stress rates, defined as the material derivative of the stress tensor minus the commutator of the spin W with stress. In the Huang VUMAT, objectivity is achieved through the incremental rotation tensor DROT: before each constitutive integration, stress is rotated from the global coordinate system to a spin-free coordinate system; after integration, it is rotated back. This method is equivalent to stress update using the Jaumann rate.

S.3 各向异性弹性响应S.3 Anisotropic Elastic Response

单晶材料的弹性响应本质上是各向异性的，其弹性刚度矩阵在晶体坐标系中通常只有少量独立常数。对于立方晶体，只有C11、C12和C44三个独立常数；对于六方晶体，有C11、C12、C13、C33和C44五个独立常数。在全局坐标系中，弹性刚度矩阵通过四阶张量变换法则进行旋转。The elastic response of single-crystal materials is inherently anisotropic, with the elastic stiffness matrix in the crystal coordinate system typically having only a few independent constants. For cubic crystals, there are only three independent constants: C11, C12, and C44. For hexagonal crystals, there are five: C11, C12, C13, C33, and C44. In the global coordinate system, the elastic stiffness matrix is rotated according to the fourth-order tensor transformation law.

在晶体塑性计算中，弹性刚度矩阵通常在子程序初始化时计算并存储，以避免每个增量步重复计算。对于小变形假设，弹性刚度矩阵可以视为常数；但对于大变形，需要考虑几何非线性对应力-应变关系的影响。In crystal plasticity computations, the elastic stiffness matrix is typically computed and stored during subroutine initialization to avoid repeated computation each increment. Under small deformation assumptions, the elastic stiffness matrix can be treated as constant; for large deformations, geometric nonlinear effects on the stress-strain relation must be considered.

S.4 状态变量持久化S.4 State Variable Persistence

在显式动力学分析中，Abaqus/Explicit会在每个增量步开始时将STATEOLD复制到STATENEW，然后调用VUMAT。VUMAT负责更新STATENEW中的值。在增量步结束时，Abaqus会将STATENEW保存为下一增量步的STATEOLD。这种设计使得状态变量的时间推进非常明确。In explicit dynamic analysis, Abaqus/Explicit copies STATEOLD to STATENEW at the beginning of each increment and then calls VUMAT. VUMAT is responsible for updating values in STATENEW. At the end of the increment, Abaqus saves STATENEW as the next increment’s STATEOLD. This design makes the time advancement of state variables very explicit.

用户需要确保在VUMAT中正确更新所有需要持久化的状态变量。任何遗漏的更新都会导致下一增量步使用错误的旧状态，从而引入累积误差。建议在子程序末尾添加状态变量完整性检查，确保所有必需变量都已更新。Users must ensure all state variables requiring persistence are correctly updated in VUMAT. Any omitted update will cause the next increment to use incorrect old states, introducing cumulative errors. It is recommended to add state variable completeness checks at the end of the subroutine to ensure all required variables have been updated.

附录G 学习路径建议Appendix G Suggested Learning Path

G.1 初学者路径G.1 Beginner Path

对于刚接触晶体塑性有限元方法的读者，建议按照以下顺序学习：（1）掌握连续介质力学基础，特别是应力、应变和本构关系的基本概念；（2）学习晶体学和滑移系理论，理解Schmid定律；（3）熟悉Abaqus用户子程序的基本架构，先掌握UMAT再过渡到VUMAT；（4）阅读黄永刚原始报告，理解代码结构；（5）运行简单的单晶拉伸算例，验证结果。For readers new to crystal plasticity finite element methods, the following learning sequence is recommended: (1) master fundamentals of continuum mechanics, especially basic concepts of stress, strain, and constitutive relations; (2) learn crystallography and slip system theory, understanding Schmid’s law; (3) become familiar with the basic architecture of Abaqus user subroutines, mastering UMAT before transitioning to VUMAT; (4) read Huang’s original report to understand code structure; (5) run simple single-crystal tension examples and verify results.

G.2 进阶路径G.2 Advanced Path

对于已经掌握基础知识的读者，可以深入研究以下内容：（1）不同硬化模型（Voce、Kocks-Mecking、Dislocation Density）的数学形式和物理背景；（2）多晶均匀化方法（Taylor、Sachs、Self-Consistent、FE2）；（3）晶体塑性中的损伤和断裂模型；（4）GPU加速和并行计算技术在晶体塑性中的应用；（5）机器学习方法在晶体塑性参数识别中的应用。For readers who have mastered the basics, the following advanced topics can be studied: (1) mathematical forms and physical backgrounds of different hardening models (Voce, Kocks-Mecking, Dislocation Density); (2) polycrystal homogenization methods (Taylor, Sachs, Self-Consistent, FE2); (3) damage and fracture models in crystal plasticity; (4) GPU acceleration and parallel computing techniques in crystal plasticity; (5) machine learning methods for crystal plasticity parameter identification.

G.3 实践项目建议G.3 Practice Project Suggestions

• 项目1：实现单晶铝的拉伸模拟，对比不同取向的应力-应变响应。Project 1: Implement single-crystal aluminum tension simulation, compare stress-strain responses for different orientations.
• 项目2：构建多晶铜的代表性体积元，模拟单轴压缩下的织构演化。Project 2: Build a polycrystalline copper representative volume element, simulate texture evolution under uniaxial compression.
• 项目3：在VUMAT中加入绝热加热效应，模拟高速剪切变形中的温度升高。Project 3: Add adiabatic heating effects to VUMAT, simulate temperature rise during high-speed shear deformation.
• 项目4：对比VUMAT和UMAT在准静态问题中的计算效率和结果差异。Project 4: Compare computational efficiency and result differences between VUMAT and UMAT in quasi-static problems.