A mechanics-informed deep learning constitutive model for sequential prediction of strain rate-dependent behavior and microstructural evolution

Gated recurrent unit

GRUs, as efficient RNN variants, have proven effective for constitutive modeling in CP simulations. Their gating mechanisms effectively address gradient vanishing issues in traditional RNNs while preserving long-term dependencies with computational efficiency. For input sequence x = (x₁, x₂, …, x_T), the GRU cell regulates information flow through update and reset gates:

where z_t and r_t represent update and reset gates, h_t denotes the hidden state, x_t is the current input, and $$ \odot $$ indicates the Hadamard product. The update gate controls previous state retention, while the reset gate modulates historical information’s influence on the current candidate state. Here, $$ W_{z}, W_{r}, W_{h} \in R^{d_{h} \times\left(d_{h}+d_{x}\right)} $$ are learnable weight matrices, $$ b_{z}, b_{r}, b_{h} \in R^{d_{h}} $$ are bias vectors, d_h denotes hidden dimension, d_x denotes input dimension, and $$ \tilde{h}_{t} $$ represents the candidate hidden state.

In classical constitutive formulations, loading paths induce discrete updates of microstructure-dependent internal variables. GRUs function as learned, discrete-time state-space models of this evolution, combining sequential interpretability with intrinsic state-update operators. By capturing temporal dependencies, they recover mechanical states governed by both instantaneous loading and accumulated deformation history. The selective retention of critical temporal features through gating mechanisms enables accurate predictions across diverse loading trajectories. Compared to traditional RNNs, GRUs achieve comparable accuracy with reduced computational cost, making them particularly suitable for the present study.

Multi-head attention mechanism

MHA, a cornerstone of modern DL architectures^[21], captures cooperative interactions among deformation mechanisms through parallel processing across multiple attention heads. Each head focuses on distinct aspects of the input sequence, extracting mechanism-specific representations across multiple timescales and loading conditions. For input sequence $$ X \in \mathbb{R}^{T \times d} $$, where T denotes timesteps and d represents mechanical feature dimensions, single-head self-attention is formulated as:

where $$ Q=X W^{Q} \in R^{T \times d_{k}} $$, $$ K=X W^{K} \in R^{T \times d_{k}} $$, $$ V=X W^{V} \in R^{T \times d_{V}} $$ are query, key, and value matrices, obtained through learnable projection matrices $$ W^{Q}, W^{K} \in R^{d \times d_{k}} $$ and $$ W^{V} \in R^{d \times d_{V}} $$, respectively. The scaling factor $$ \sqrt{d_{k}} $$ controls dot product magnitudes to prevent gradient vanishing.

To enhance expressivity and capture inherent features tied to dislocation and twinning, MHA projects inputs into h subspaces, computes self-attention in each, then concatenates results with linear transformation:

where head_i = Attention (Q_i, K_i, V_i), W_o denotes the output transformation matrix, and h represents the number of attention heads. The input feature dimension must be divisible by h; if not, a 1D convolution can adjust it. Individual heads may specialize in specific mechanisms (dislocation dynamics, strain-rate sensitivity, or coupled effects), enabling simultaneous focus on distinct representation subspaces and enhancing feature diversity and multiscale characterization.

Residual connections enhance information flow and mitigate gradient vanishing:

Layer normalization stabilizes training through feature normalization:

where d represents feature dimensions, $$ \epsilon_{s} $$ ensures numerical stability, μ_LN and σ_LN denote the mean and standard deviation of the input features used in the layer-normalization operation, $$ \epsilon_{s} $$ is a small constant for numerical stability, and γ, β are learnable scaling and bias parameters. The normalization suppresses internal covariate shifts and accelerates convergence. The integration of MHA, residual connections, and layer normalization constructs a stable and efficient network architecture.

Mechanics-informed layer

The proposed mechanics-informed layer integrates classical elasto-plastic theory with deep learning to regularize predicted stress-strain responses. Derived from conventional J2 flow theory, this layer is implemented through differentiable neural network formulations that preserve end-to-end differentiability while incorporating fundamental mechanics principles. By assuming initially random crystallographic texture, it enables differentiable mapping between principal strains and stresses, with tensors expressed in principal form.

Since crystal plasticity involves continuous texture evolution and lattice rotation, fixed-axis anisotropic criteria (e.g., Hill48) are unsuitable. Instead, the evolving anisotropy is captured by GRU-predicted residual terms, effectively enabling the model to learn a dynamic yield surface. It should be noted that the J₂ plasticity formulation within the mechanics-informed layer serves strictly as a universal isotropic elastoplastic baseline, rather than a complete description of evolving anisotropy. This architectural decoupling leverages the principle of residual learning: the physics layer governs the fundamental isotropic transition, freeing the data-driven components to fully dedicate their representational capacity to extracting higher-order anisotropic residuals. The effectiveness of this approach in breaking the isotropic assumption is explicitly evidenced by the model’s ability to capture the non-linear distortion of the yield surface under varying multiaxial loading ratios, demonstrating its capability in representing pronounced directional dependencies. While the J₂ baseline inherently assumes tension-compression symmetry, the data-driven residual terms effectively compensate for this limitation, allowing the model to capture the distinct tension-compression asymmetry induced by the polarity of deformation twinning at high strain rates.

The layer processes a 3D principal strain tensor ε = [ε₁, ε₂, ε₃] as input. The initial estimate of the corresponding stress tensor, σ_e,i is computed using a simplified isotropic elasticity model:

where λ and μ are Lamé constants, treated as learnable parameters within the model. To evaluate potential plastic deformation, the von Mises yield criterion is implemented. The equivalent trial stress, $$ \bar{\sigma}_{e} $$, is calculated as:

where $$ \sigma_{e, i}^{\prime}=\sigma_{e, i}-\frac{1}{3}\left(\sigma_{e, 1}+\sigma_{e, 2}+\sigma_{e, 3}\right) $$ is the deviatoric tensor component of the trial stress. The yield function ϕ is defined as $$ \phi=\bar{\sigma}_{e}-\sigma_{y} $$, where σ_y is implemented as a learnable parameter. To ensure differentiability across the elastic-plastic transition, a smooth transition is achieved by using a sigmoid function:

with parameter k controlling the transition steepness, effectively regularizing the yield surface sharpness.

During plastic deformation (ϕ > 0), the plastic strain increment Δλ is determined as:

where H is a learnable hardening parameter accounting for strength evolution due to plastic strain. The corresponding stress is subsequently updated to include plastic contributions:

Consequently, the layer outputs three-dimensional (3D) principal stress tensor σ = [σ₁, σ₂, σ₃]. This mechanics-informed approach preserves classical plasticity foundations while providing adaptive capabilities through physics-constrained learnable parameters, simultaneously offering physical interpretability and data-driven flexibility for material behavior simulation.

Mechanics-informed loss function

To address the mechanics-informed layer requirements, a mechanics-informed loss function is developed, designed to ensure both data fitting accuracy and physical fidelity. The multi-term loss function is formulated as:

The foundational mean squared error (MSE) loss quantifies the deviation between predicted and ground truth values:

where N denotes the total sample count, and y_i,pred, y_i,true represent predicted and actual mechanical quantities (stress components, DD, TVFs, ensuring accuracy across strain-rate dependent responses and microstructural features.

To capture temporal evolution with physical continuity, a temporal consistency constraint loss $$ \mathcal{L}_{\text {temporal }} $$ is implemented:

where y_t and y_t+1 represent consecutive time-step quantities. This term penalizes discontinuities in predicted mechanical quantity evolution, enforcing continuous deformation laws.

For optimized performance in the elastic regime, prevalent in the initial stages of crystal plasticity simulations, an elastic region constraint loss is implemented. This specialized term ensures that the model exhibits linear elastic behavior under small strain conditions, adhering to Hooke’s law:

where $$ \mathcal{L}_{\text {temporal }} $$ is an indicator function equaling 1 within the elastic regime, and 0 otherwise. The term 12 prevents numerical instability. Here, the elastic constraint is explicitly defined using a linear stress-strain proportionality to enforce Hooke-type behavior in the small-strain regime.

The yield criterion loss enforces physically consistent elastic-plastic transition behavior, constraining stress predictions to remain within physically meaningful boundaries:

By incorporating the yield stress σ_y, hardening coefficients H, and plastic strain ε_p,i, the model accurately simulates the nonlinear mechanical responses. This component ensures that predicted stress trajectories respect fundamental plasticity principles, preventing unphysical extrapolations.

The comprehensive loss function is formulated as a weighted sum of individual components, with coefficients λ₁ - λ₄ adjustable to specific applications and optimization objectives. This parameterization allows precise modulation of model focus, increasing λ₃ enhances elastic regime performance, while amplifying λ₂ strengthens physical plausibility constraints.

These embedded mechanical constraints provide intrinsic regularization, reducing training data requirements, which is critical for data-scarce materials research. The MIDLCM transforms traditional data-intensive paradigms into knowledge-enhanced learning processes, leveraging established physical principles for accurate, physically consistent, and empirically effective regularization. Note that in the present study, the specific weight allocation (λ₁ = 1.0, λ_2-4 = 0.1) was determined based on physical intuition rather than an exhaustive sensitivity analysis. This 10:1 ratio establishes a hierarchical regularization strategy: the data-fitting term (λ₁) acts as the primary optimization objective, while the mechanics-informed terms (λ_2-4) serve as auxiliary soft constraints. As validated by the subsequent ablation study, this empirically established ratio effectively enhances physical consistency without overwhelming the data-driven learning process, representing a reasonable baseline setting rather than a mathematically optimal choice. These current weights establish a stable hierarchical relationship; nevertheless, precisely optimized values may be further determined using adaptive loss weighting techniques.

A brief introduction of training dataset generator (EVPSC model)

EVPSC model^[30] is employed for its comprehensive integration of multiple deformation mechanisms with strain-rate sensitivity while maintaining computational efficiency superior to conventional crystal plasticity models, enabling practical training dataset generation within reasonable computational costs. The EVPSC model incorporates dislocation and twinning mechanisms within a self-consistent crystal plasticity formulation^[1,31]. The constitutive formulation describes the elastic-plastic deformation of individual grains through the deformation gradient tensor, where elasticity is characterized by the Jaumann rate of Cauchy stress coupled with anisotropic stiffness tensors^[30].

We incorporated multiple mechanisms for a reliable description of the strain rate sensitive mechanisms and mechanical behaviors. Dislocation-mediated plasticity follows Orowan’s equation^[32], with dislocation dynamics (nucleation, multiplication, annihilation) governed by differential equations contributing to strain-rate sensitivity^[33,34]. Concurrently, the TDT model captures twin nucleation and evolution^[30]. Grain boundaries influence plasticity through the Hall-Petch effect^[35-39]. Dislocation pile-up and boundary nucleation are incorporated into the DD equations, dynamically updating slip resistances and enabling strain-rate-dependent hardening^[40,41]. The self-consistent framework homogenizes grain interactions by iteratively solving stress equilibrium and compatibility across the polycrystal, simulating uniaxial tension at different strain rates^[31]. This approach captures macroscopic stress-strain responses, texture evolution, and TVF changes, with time-stepping algorithms resolving the coupling between grain-scale mechanisms and polycrystalline homogenization.

The polycrystalline system comprises 5,000 grains with 90% face-centered cubic (FCC) matrix strengthened by 10% body-centered cubic (BCC) particles. BCC particles are incorporated by increasing slip resistances on 12 {111}<110> slip systems in the FCC phase, approximating their barrier effect on dislocation motion without explicit mesoscale modeling. Strain rates range from 0.001/s to 4,000/s. {111}<112> twinning also contributes to plasticity, becoming more significant at high strain rates, consistent with experimental observations. A more detailed formulation, the calibrated model parameters, and the nomenclature are provided in the Supplementary Materials [Supplementary Tables 1 and 2]. The EVPSC dataset assumes linear elasticity, justifying the Hooke’s Law constraint. For extensions into shock physics, this constraint requires modification to accommodate adiabatic heating and shock-induced nonlinear elastic effects.

Parameter calibration of EVPSC for CrFeNi alloys

To establish a reliable dataset for training the constitutive model, numerical simulations are performed using the aforementioned EVPSC model, with subsequent validation against experimental CrFeNi alloy observations. Model parameters are calibrated using tensile deformation data across multiple strain rates, with experimental validation encompassing four strain rates (0.001/s, 0.1/s, 3,000/s, and 4,000/s) along the rolling direction^[26]. The initial microstructure features a random crystallographic texture (see Figure 4A). Parameter optimization yields an optimal constitutive parameter set (see Supplementary Table 1), enabling high-fidelity dataset generation. The calibrated model accurately reproduces CrFeNi’s mechanical response across the deformation regime, with predictions showing excellent agreement with experimental stress-strain data [Figure 4B], as well as the corresponding DD evolution [Figure 4C] and TVF evolution [Figure 4D] at different strain rates. This consistency reflects the model’s accurate representation of fundamental mechanisms, including dislocation nucleation, multiplication, grain boundary accumulation, and strain hardening.

Figure 4. (A) Electron backscatter diffraction map illustrating the initial texture of the CrFeNi alloy prior to deformation, depicted through an inverse pole figure color scheme^[26]; (B-D) Comparison of the mechanical response predictions including stress, DD and twin volume fraction derived from crystal plasticity (CP) simulations to experimental (Exp.) measurements obtained at varying strain rates. Experimental results adapted from Wang et al.’s study^[26].

The EVPSC predictions exhibit clear strain-rate-sensitive behavior. The strain rate sensitivity of the yielding stress is markedly higher near 3,000/s than under quasi-static conditions, and only the dynamic cases show a rapid initial rise in DD and activation of twinning. Accurately reproducing these strain rate sensitive features is essential yet challenging for deep learning constitutive models, and our proposed MIDLCM is designed to capture them. This reproduction will underscore the advantages of the present model design.

High-fidelity dataset generation and pre-processing

High-fidelity training data are generated using the calibrated CrFeNi polycrystalline model under plane stress conditions, with specimens subjected to biaxial tension at loading angle $$ \alpha=\arctan \left(\frac{\varepsilon_{2}}{\varepsilon_{1}}\right) $$, as shown in Figure 3. Remaining strain components satisfy stress-free boundary conditions. As shown in Figure 4A, a single microstructure with random crystallographic orientations is employed across multiple loading directions to optimize computational efficiency while ensuring dataset representativeness. Simulations span six strain rates (0.001/s, 0.1/s, 10/s, 1,000/s, 3,000/s, and 5,000/s) covering quasi-static to high rate regimes, with maximum equivalent strain limited to 0.16. Critical micromechanical features, including stress-strain tensors, dislocation densities, and TVFs, are recorded at each increment. The dataset employs a single random-texture microstructure, with 756 cases spanning six strain rates, and three crystallographic planes. Rather than an algorithmic limitation, this approach serves as a deliberate scientific control. The primary objective of this framework is to isolate and project the complex rate-dependent kinetics - specifically, the transition from dislocation-dominated plasticity to twinning-assisted mechanisms under extreme dynamic loading - into the deep learning architecture. Simultaneously introducing variations in initial grain size, texture, and phase fractions would heavily convolute these geometric and crystallographic variables with the kinetic rate-dependent mechanisms, thereby obscuring the underlying physical causality. By maintaining a controlled initial microstructural state, we ensure that the framework’s capability to accurately capture explicit macro-micro coupled physics across varying strain rates and loading paths can be rigorously and unambiguously evaluated.

Traditional constitutive models employ sequential incremental integration, computing stress and microstructural states step-by-step through iterative time-stepping procedures, as shown in Figure 1B. Conversely, DL constitutive models adopt parallel computation, accepting complete strain histories as input and simultaneously predicting full stress responses and microstructural evolution in a single forward pass, achieving orders of magnitude computational acceleration while eliminating incremental procedures.

Data preprocessing enhances model convergence and performance. Logarithmic transformation is applied to dislocation densities and strain rates to compress multi-order magnitude variations and improve model sensitivity across scales. Subsequently, all features, including strain-stress histories, strain rates, dislocation densities, and TVFs, undergo standardization to zero mean and unit variance:

where μ and σ represent the mean and standard deviation of the original data.

The dataset is partitioned into training and testing sets at a 4:1 ratio, with representative sampling ensuring balanced distributions across strain paths and rates. The case-based splitting assigns each case, defined by strain rate, loading angle, and crystallographic plane, atomically to training or test sets, preventing temporal correlation leakage while ensuring test cases represent genuinely unseen loading paths.

Model training configuration

All experiments are conducted on a computer equipped with an NVIDIA RTX 3090 GPU, utilizing PyTorch 2.0.0 and CUDA 11.8 to leverage GPU-accelerated computation. The optimization strategy is meticulously designed to ensure model convergence, stability, and performance. The AdamW optimizer^[42] is employed, combined with a cosine annealing learning rate scheduler. AdamW, an advanced variant of the Adam, decouples weight decay from gradient updates, enhancing generalization through explicit weight regularization.

The hyperparameters for AdamW are configured with β₁ = 0.9, β₂ = 0.999, a weight decay coefficient of 0.01, an epsilon value of 10^-8, and a learning rate of 0.001. A cosine annealing scheduling dynamically adjusted the learning rate throughout training. This scheduler gradually increases the learning rate during initial training to prevent parameter instability, subsequently modulating the learning rate through a cosine annealing function:

where η_min and η_max represent minimum and maximum learning rates, t denotes the current step, and T represents the total steps. The first 10% of steps are allocated for warm-up, with η_max = 0.001, η_min = 0, and T = 1,000. Compared to traditional linear decay methods, this scheduler enables effective parameter space exploration while mitigating convergence to suboptimal local minima.

Model weights are initialized using Xavier initialization to establish appropriate initial parameter distributions. MIDLCM is trained from scratch for 1,000 epochs with a batch size of 64. Gradient clipping with a maximum norm threshold of 1.0 is implemented to enhance training stability and prevent gradient explosion. Validation performance is continuously monitored during training, with optimal model parameters preserved based on validation loss minimization. This protocol ensured robust convergence while maintaining physical consistency in the learned representations of material behavior. Training stability was promoted through Xavier initialization, cosine-annealed AdamW optimization, gradient clipping, and validation-based checkpoint selection.

The present MIDLCM is formulated as a deterministic constitutive framework, and no formal uncertainty quantification is considered in this study. Variance-based sensitivity analysis and probabilistic treatment of uncertain material parameters or loading conditions remain important directions for future work when extending the model to reliability-oriented engineering applications.

Hyperparameter sensitivity analysis

This study employs Bayesian optimization with Tree-structured Parzen Estimator^[43] for efficient hyperparameter tuning of MIDLCM. Unlike traditional grid search methods, Bayesian optimization constructs a probabilistic model of the objective function (validation MSE), and enables directed exploration of the high-dimensional parameter space through sequential sampling with acquisition functions that balance exploitation and exploration. As shown in Figure 2, the hyperparameter search space encompassed: GRU hidden dimensions, GRU layer number N_G ∈ {1, 2, 3, 4}, embedding dimensions d_E ∈ {1, 2, 3, 4}, MHA heads h ∈ {1, 2, 4, 8}, and FFN hidden dimensions d_F ∈ {64, 128, 256, 512}. Optimization converged after 35 iterations out of 50 total iterations, yielding an optimal configuration of a two-layer GRU with 256 hidden dimensions, three embedding dimensions, four attention heads, and an FFN hidden dimension of 128. This configuration balances representational capacity with computational efficiency while minimizing overfitting risk.

As shown in Figure 5, sensitivity analysis results reveal distinct performance patterns across architectural configurations. GRU hidden dimensions exhibit a non-monotonic relationship with prediction accuracy, with MSE decreasing from 0.0386 to 0.0282 as dimensions increase from 64 to 256, then deteriorating at 512 dimensions (MSE = 0.0341), indicating an optimal representational capacity threshold beyond which overfitting occurs. Network depth analysis confirms optimal performance at two GRU layers, with both shallower and deeper architectures showing reduced efficacy. Embedding dimensionality achieves peak performance at three dimensions, sufficiently capturing the rate feature space complexity. MHA mechanism peaks at four heads, indicating that moderate attention complexity best facilitates feature integration. FFN performance maximizes at 128 hidden dimensions, with larger configurations degrading due to overparameterization.

Figure 5. Experimental results of hyperparameter sensitivity. (A-E) MSE variations across key hyperparameters: GRU hidden dimensions, GRU layers, embedding dimensions, MHA heads, and FFN hidden dimensions; (F) Normalized importance scores of investigated hyperparameters. MSE: Mean squared error; GRU: gated recurrent unit; MHA: multi-head attention; FFN: feed-forward network.

These findings demonstrate that model performance maximizes at intermediate complexity levels across all hyperparameters, emphasizing the critical balance between representational capacity and overfitting in crystal plasticity constitutive modeling. The narrow MSE variation range (0.028-0.041) indicates robust architectural performance with limited hyperparameter sensitivity, suggesting inherent stability in capturing crystal plasticity mechanisms. Normalized importance analysis reveals that for the rate-dependent modeling task, embedding dimension exhibits the highest relative importance, followed sequentially by MHA heads, FFN hidden dimensions, GRU hidden dimensions, and GRU layer number.

Ablation experiments

Ablation experiments of MIDLCM quantify individual component contributions by progressively incorporating MHA, mechanics-informed layer, and mechanics-informed loss function into a baseline GRU network. Performance metrics include MSE, mean absolute error (MAE), parameter count, and inference time, with results summarized in Figure 6.

Figure 6. Performance and computational cost results of model ablation experiments examining the four principal components of MIDLCM architecture. MSE: Mean squared error; MAE: mean absolute error; GRU: gated recurrent unit; MHA: multi-head attention; MIDLCM: mechanics-informed deep learning constitutive model.

The baseline GRU network achieved moderate accuracy (MSE = 0.062, MAE = 0.184) with computational efficiency (5.17 ms inference, 0.74M parameters). Incorporating MHA substantially improved performance, reducing MSE and MAE by 42.0% and 37.5% respectively, with minimal computational overhead (1.00M parameters, 5.75 ms). This improvement demonstrates MHA’s effectiveness in capturing temporal dependencies inherent in mechanical and microstructural evolution. By mitigating the vanishing-gradient issues of standard GRUs, MHA ensures that critical historical shifts across the extreme strain-rate spectrum are accurately preserved and transmitted to the subsequent mechanics-informed layer, laying the mathematical foundation to enforce physical consistency. Sequential integration of the mechanics-informed layer and loss function yielded progressive enhancements, ultimately achieving optimal performance (MSE = 0.028, MAE = 0.100) while maintaining computational efficiency (5.91 ms). The mechanics-informed components contributed incremental but consistent improvements of 13.9% and 9.7% in MSE reduction, respectively. These incremental gains confirm that the improvement arises from the powerful synergy of explicitly tracking rate-dependent mechanism shifts, mechanics-informed elastoplastic regularization, and sequence interaction learning, rather than from using the GRU architecture alone. To further evaluate algorithmic robustness against data partitioning and local minima, comprehensive five-fold cross-validation and multi-seed repeated runs were conducted on the optimized MIDLCM framework. The five-fold case-based cross-validation yielded a mean test MSE of 0.0288 ± 0.0010 and a mean test MAE of 0.1016 ± 0.0024, indicating stable generalization under different train-test splits. Furthermore, five independent training runs using different random initialization seeds resulted in a mean MSE of 0.0282 ± 0.0005 and MAE of 0.1002 ± 0.0013. These results, detailed in Supplementary Materials [Supplementary Table 3], explicitly demonstrate that the framework is highly stable, maintaining consistent predictive performance without acute sensitivity to random stochasticity.

These results demonstrate the synergistic effect of combining deep learning with mechanics-informed constraints. Each component addresses distinct aspects of crystal plasticity modeling: MHA captures complex temporal patterns, while mechanics-informed components embed domain knowledge into both architecture and optimization. The progressive improvement indicates that the mechanics-informed components improve predictive accuracy and physical consistency relative to the purely data-driven baseline for constitutive modeling, achieving superior accuracy without computational penalty through principled integration of mechanical principles.

Computational efficiency analysis

Evaluating computational efficiency requires careful context. While EVPSC serves as the high-fidelity ground truth resolving comprehensive grain-level details, MIDLCM significantly accelerates computations by capturing the overall microstructural evolution. For single-case predictions, MIDLCM processes 100-timestep strain paths in 0.005 s, while EVPSC requires approximately 120 s - a 24,000-fold speed-up. This computational advantage is further enhanced in batch processing. For parallel predictions across 100 distinct loading conditions, MIDLCM completed analyses in 0.07 s compared to approximately 200 min for EVPSC. Furthermore, comparing MIDLCM directly to purely macroscopic deep learning models would be unequal, as it overlooks the inherent computational complexity of explicitly tracking internal state variables like DD and TVF. Ultimately, this acceleration preserves critical mechanism-level fidelity, establishing the framework as highly efficient for rapid material behavior assessment and multi-scale simulations.

While empirical macroscopic models (e.g., Johnson-Cook) offer extreme computational efficiency for industrial FEM simulations, they inherently trade microstructural fidelity for speed and cannot spontaneously capture dynamic state-driven mechanism transitions. Achieving comparable microstructural fidelity requires traditional constitutive frameworks with multiple internal state variables and implicit integration of highly nonlinear, stiff evolution equations (ODEs), substantially increasing computational costs. In contrast, MIDLCM avoids these expensive iterative updates through a unified deep learning architecture, providing an efficient surrogate while retaining explicit microstructural descriptors. Because the present study focuses on mechanism-coupled predictions rather than purely macroscopic stress mapping, a direct benchmark against calibrated traditional phenomenological models is not included. Therefore, the current computational results should be strictly interpreted as a highly favorable fidelity-efficiency tradeoff relative to the EVPSC reference model, rather than a general superiority claim over all macroscopic constitutive formulations.

Uniaxial tension prediction competence compared to experimental data

The proposed MIDLCM predictions for uniaxial tension at strain rates of 0.001/s, 0.1/s, and 3,000/s are validated against experimental data from Wang et al.^[26]. As shown in Figure 7A, MIDLCM demonstrates excellent agreement with the experimental stress-strain curves across the examined strain rates. The model accurately captures microstructural evolution, including initial dislocation proliferation during early plasticity and subsequent saturation due to dynamic recovery processes, as shown in Figure 7B. TVF predictions effectively reproduce strain rate-dependent activation and growth patterns, demonstrating the model’s capability to represent underlying deformation mechanisms across multiple orders of magnitude in strain rate.

Figure 7. Evolution of stress, DD, and TVF with deformation in polycrystalline CrFeNi samples with strain rates of 0.001/s, 0.1/s, and 3,000/s. TVF: Twin volume fraction; DD: dislocation density; DL: deep learning.

To evaluate extrapolation performance, MIDLCM is tested on CrFeNi alloys under uniaxial tension at strain rates of 0.0001/s, 2,000/s, and 4,000/s - conditions outside the training dataset. Figure 8 illustrates remarkable consistency between predicted stress-strain curves and experimental measurements across these untrained conditions. The model accurately captures strain rate-dependent behavior, including elastic response, yield phenomena, and strain hardening characteristics. Microstructural predictions successfully reproduce DD and TVF evolution from quasi-static (0.0001/s) to ultra-high strain rate conditions (4,000/s), demonstrating robust representation of rate-dependent mechanisms beyond training boundaries. This extrapolation capability validates the physical consistency of the mechanics-informed architecture, indicating a genuine understanding of underlying physics rather than data interpolation. The demonstrated accuracy and computational efficiency across diverse loading conditions establish the model as an effective alternative to computationally intensive crystal plasticity simulations, particularly valuable for engineering applications with limited experimental data availability.

Figure 8. Evolution of stress, DD, and TVF with deformation in polycrystalline CrFeNi samples with strain rates of 0.0001/s, 2,000/s, and 4,000/s. TVF: Twin volume fraction; DD: dislocation density; DL: deep learning.

Prediction competence under biaxial loading conditions

Quantitative validation demonstrates excellent agreement between MIDLCM predictions and crystal plasticity simulations. Figure 9 presents parity plots for timestep-averaged equivalent stress, dislocation density, and TVF under biaxial loading conditions, where predicted values closely align with reference solutions along the y = x diagonal across all samples. This strong correlation validates the model’s capability to simultaneously capture macroscopic mechanical response and underlying microstructural evolution with high fidelity.

Figure 9. Parity plots comparing MIDLCM predictions with reference dataset under biaxial loading conditions for (A) time-averaged equivalent stress; (B) time-averaged DD; and (C) time-averaged twin volume fraction. MIDLCM: Mechanics-informed deep learning constitutive model; DD: dislocation density.

To demonstrate MIDLCM’s predictive accuracy across varying performance levels, three representative samples with different MSE magnitudes and loading angle α across different planes are examined in Figure 10. These samples are strategically selected to showcase predictions spanning from high-accuracy cases to moderate-precision scenarios and more challenging predictions, collectively illustrating the model’s capability range across diverse loading conditions. Compared to reference dataset, MIDLCM accurately reproduces elastic behavior, yield transition, and strain hardening characteristics. Low strain rates exhibit gradually declining hardening rates, while high rates induce initial stress fluctuations followed by rapid stabilization and elevated hardening. These variations reflect strain rate-dependent activation of dislocation and twinning mechanisms.

Figure 10. MIDLCM (labeled as DL) and reference dataset (labeled as Ref.) results for the mechanical response of CrFeNi alloys with strain rates from 0.001/s to 5,000/s. MIDLCM: Mechanics-informed deep learning constitutive model; DL: deep learning.

DD evolution demonstrates the model’s microscale accuracy. The model captures nonlinear behavior, including initial rapid growth and subsequent saturation. Below 2,000/s, DD increases steadily; above this threshold, accelerated proliferation occurs at plastic flow onset. This behavior, governed by grain boundary nucleation mechanisms, is essential for representing dislocation hardening under high strain rates. Higher saturation densities at elevated rates result in increased final dislocation densities.

For twinning behavior, MIDLCM correctly predicts negligible twin nucleation under quasi-static loading and significant activation at medium-to-high strain rates. Twin proliferation rates vary with strain rate due to stress-dependent growth kinetics: below 10/s, twinning propagation remains minimal, while higher rates induce substantial acceleration. MIDLCM successfully reproduces these strain rate-dependent characteristics, capturing the complex interplay between coupled deformation mechanisms and validating its capability for strain rate-sensitive phenomena simulation.

Twinning mechanisms significantly influence mechanical behavior through orientation-dependent deformation that induces asymmetric responses. To investigate this phenomenon, yield surfaces are predicted under biaxial loading at 0.001/s and 3,000/s across equivalent plastic strains of 0.01, 0.05, 0.1, and 0.15. Data are extracted from xy, yz, and xz principal strain planes in the sample coordinate system (x: rolling direction, z: normal direction) as shown in Figure 11.

Figure 11. Prediction of yield surface shape with different plastic strains for CrFeNi polycrystalline system with strain rates of 0.001/s and 3,000/s. 2: strain rate.

Under quasi-static conditions (0.001/s), yield surfaces maintain symmetry throughout deformation, exhibiting a uniform expansion characteristic of isotropic hardening dominated by dislocation-mediated plasticity. At dynamic rates (3,000/s), initial yield surfaces remain symmetric, but subsequent deformation induces pronounced asymmetry due to twinning activation. This demonstrates how the data-driven residual terms compensate for the symmetry of the J₂ baseline to capture the polar nature of dynamic twinning. Comparative analysis reveals distinct morphological differences and substantial stress magnitude variations between strain rates. These results demonstrate strain rate sensitivity in various aspects, encompassing concurrent changes in flow stress, yield surface morphology, mechanical anisotropy, and tension-compression asymmetry - all accurately captured by MIDLCM. Due to the extreme challenges of dynamic multiaxial experiments, these multiaxial predictions currently constitute a simulation-to-simulation benchmark. However, the rigorously calibrated EVPSC model provides a physically reasonable reference, with its true accuracy awaiting future experimental confirmation.

Furthermore, MIDLCM predictions are compared with reference dataset to assess dislocation density and TVF evolution under varying strain rates and loading orientations [Figure 12]. Analysis is performed at equivalent strains of 0.01, 0.05, 0.1, and 0.15 to capture progressive microstructural evolution. Mechanistically, the predicted rate sensitivity originates from the dynamic competition between dislocation and twinning. Under quasi-static conditions, hardening is dominated by progressive dislocation accumulation. Conversely, at dynamic rates (3,000 s^-1), conventional dislocation mobility is severely restricted, driving a rapid, drastic defect accumulation near yielding. This physical phenomenon is consistent with dynamic mechanisms observed in highly rate-sensitive metallic systems like SAC305 solders^[44]. Explicitly linking this dislocation structure evolution to macroscopic strength ensures robust dynamic hardening predictions, aligning with recent mechanistic studies on Ni-based superalloys^[45]. Subsequently, this restricted dislocation mobility triggers extensive deformation twinning to accommodate the dynamic strain, intrinsically governing the non-linear evolution and pronounced asymmetry of the yield surfaces (as observed in Figure 11). Capturing such profound mechanism mutations across critical rate thresholds is universally essential for dynamically loaded materials, a necessity similarly highlighted in extreme thermo-mechanical studies of biomimetic polymers^[46]. Ultimately, this quantitative agreement validates the framework’s reliability in capturing complex dislocation-twinning interplays.

Figure 12. Prediction of dislocation density and TVF with different plastic strains for the CrFeNi polycrystalline system under varying strain rates and load conditions. CP: Crystal plasticity; DL: deep learning; TVF: twin volume fraction.

Visualization analysis of MHA under different strain rates

Deep learning models often encode mechanisms in ways that are hard to interpret. As shown previously, MHA substantially improves predictive performance by tracking temporal trends, visualized via average attention heatmaps and diagonal evolution curves in Figure 13. MIDLCM’s interpretability fundamentally derives from its mechanics-informed layer enforcing plasticity principles, alongside explicit predictions of microstructural descriptors (dislocation density, TVF). Therefore, these MHA visualizations serve strictly as diagnostic tools revealing the network’s temporal focus across loading stages.

Figure 13. MHA-based attention visualization under different strain rates: (A) attention heatmaps; and (B) diagonal attention evolution. MHA: Multi-head attention.

In the heatmaps, color intensity denotes attention weight, with bright regions highlighting historical strain states referenced for predicting current stress. The maps are clearly strain-rate sensitive. Elevated attention around the elastic-plastic transition, particularly at lower rates, suggests yielding dictates the macroscopic response and history dependence. Moreover, early attention re-concentration mitigates GRU vanishing-gradient issues for long sequences, improving prediction fidelity. This evolving attention distribution aligns broadly with yielding and hardening behaviors. However, it must be emphasized that these learned components provide Supplementary Materials about how the model distributes attention across deformation stages; they should be interpreted only as model-level temporal attribution, rather than direct evidence of mechanism-specific causality.

An additional insight from the attention maps is that diagonal attention evolution highlights critical physical events and the model’s focus on instantaneous states during deformation. Attention variations near the yield point reflect emphasis on plastic deformation initiation, while patterns during uniform plastic deformation capture work hardening behavior. Significant attention changes at failure points validate adaptive attention adjustments during crucial physical transitions. Overall, diagonal attention values increase with strain rate, indicating prioritization of instantaneous states at high strain rates. This behavior aligns with physical intuition that high strain rate responses are more rate-dependent than history-dependent. Results demonstrate adaptive attention adjustments based on loading conditions: emphasizing instantaneous effects at high strain rates while comprehensively considering historical deformation paths at low strain rates. This attention distribution provides unique insights in MIDLCM, distinct from conventional models.

While averaged attention patterns demonstrate physically meaningful behavior, individual attention heads do not form one-to-one mappings with dislocation or twinning processes, rather, MHA primarily improves temporal modeling, and the model’s physical interpretability fundamentally derives from its explicit microstructural outputs. MHA operates as a distributed representation learning mechanism capturing complex temporal dependencies, with explicit physical constraints encoded through the mechanics-informed layer rather than pre-assigned mechanism-specific attention roles. We note that attention weights provide supplementary insight into model behavior but should be interpreted cautiously, as attention patterns do not necessarily reflect causal relationships. The primary interpretability of MIDLCM derives from its explicit microstructural outputs and mechanics-informed architecture.

The analysis of mechanics-informed components

Given that the mechanics components (mechanics-informed layer and loss) clearly enhance the model’s competence in the ablation experiments (see Figure 6), we should further emphasize their importance, particularly how they influence the outcomes. To highlight this, we built a model without the mechanics-informed components to show that this ablated model’s performance lacks detailed accuracy. The ablated model retains identical hyperparameters as the complete MIDLCM to ensure fair comparison of the mechanics-informed components’ contribution.

The results demonstrate the critical importance of mechanics-informed guidance. Compared to the complete model [Figure 9], the parity plots in Figure 14A reveal substantial performance degradation, characterized by increased scatter and deviations from the ideal diagonal. To make this degradation more intuitive, we selected a biaxial case (yz, α = 63.43°). The ablated model yields superficially plausible predictions, but the details are insufficient for quantitative, theory-driven use. Most notably, the ablated model exhibits significant deficiencies in capturing the elastic-plastic transition regime [Figure 14B], demonstrating poor accuracy at the yield point where complex dislocation mechanisms initiate. In addition, the predicted flow stress fluctuates at higher strains, which is a prevailing phenomenon across all cases. These deficiencies also impair predictions for experimental uniaxial tension over the wide range of strain rates [Figure 14C].

Figure 14. Performance of MIDLCM without mechanics-informed component: (A) The parity plot of the time-averaged equivalent stress; (B) results from a biaxial loading case (yz, α = 63.43°) conducted by the ablated model compared to the reference CP results; (C) uniaxial tension cases prediction from the ablated model compared with experiments results. MIDLCM: Mechanics-informed deep learning constitutive model; CP: crystal plasticity.

The mechanics-informed components are embedded in this framework to mitigate these deficiencies. The poor prediction near the yield point largely stems from limited elastic data, with the elastic mechanical responses available only up to 0.2% strain. To address this, the mechanics-informed layer explicitly encodes the elastic-plastic transition, and the loss includes the yield point and the elastic part. We also observe flow-stress fluctuations at high strain from the ablated model. To mitigate this, we incorporate dislocation density and twinning fraction into the final stress prediction immediately before the mechanics-informed layer. This design brings microstructural information into the stress response and yields a more stable flow-stress prediction.

As shown in Figure 10, comparing with MIDLCM indicates that the mechanics-informed components successfully overcome these deficiencies, indicating that the mechanics-informed layers effectively incorporate fundamental physical principles rather than merely overfitting training data. The mechanics-informed architecture enables accurate prediction of complex mechanical behavior while maintaining physical consistency across diverse loading conditions, demonstrating its essential role in physics-informed machine learning frameworks for materials modeling. The improvement validates that the proposed mechanics-informed neural architecture, whose principles are detailed in Sections Mechanics-informed layer and Mechanics-informed loss function, significantly outperforms purely data-driven approaches for constitutive modeling.