2.1 Johnson-Cook constitutive model

The JC constitutive model incorporates strain hardening, strain-rate effects, and thermal softening, with the flow stress defined as:$\sigma =\left(A+B{\left({\epsilon }^{\text{p}}\right)}^n\right)\left(1+C\ln \frac{{\dot{\epsilon }}^{\text{p}}}{{\dot{\epsilon }}_0}\right)\left(1-{\left(T*\right)}^m\right),$(1)where σ denotes the effective flow stress, ε^p is the effective plastic strain, ${\dot{\epsilon }}^{\text{p}}$ represents the effective plastic strain rate, and ${\dot{\epsilon }}_0$ is the reference strain rate. Parameters A, B, n, C, and m denote the initial yield stress, strain hardening coefficient, hardening exponent, strain rate sensitivity coefficient, and the thermal softening exponent, respectively. T* is the normalized temperature, which is given as:$T*=\frac{T-T_{\text{r}}}{T_{\text{m}}-T_{\text{r}}},$(2)where T denotes the current temperature, T_r denotes the room temperature, and T_m represents the melting temperature of the material. In this work, temperatures are considered to be uniform, i.e., T* = 0. Equation (1) is therefore simplified as follows:$\sigma =\left(A+B{\left({\epsilon }^{\text{p}}\right)}^n\right)\left(1+C\ln \frac{{\dot{\epsilon }}^{\text{p}}}{{\dot{\epsilon }}_0}\right).$(3)Assume that v=(v₁, …, v_N) represents an N-dimensional vector of internal state variables, the Helmholtz free energy can be defined as:$\psi =\psi \left(\mathit{\epsilon },v\right),$(4)where ψ is the Helmholtz energy potential, and ε denotes the strain vector. According to thermodynamic principles, the following relation can be derived:$\mathit{\sigma }=\frac{\partial \psi }{\partial \mathit{\epsilon }},$(5)$D=-{\displaystyle \sum _{i=1}^N\frac{\partial \psi }{\partial v_i}}\cdot {\dot{v}}_i,$(6)in which σ denotes the stress tensor and D is the dissipation rate. Both the Helmholtz energy potential and the dissipation rate are required to be non-negative in accordance with thermodynamic principles:$\psi \ge 0,$(7)$D\ge 0.$(8)

2.2 New developments

Figure 1A shows the architecture of the PIRNN model, which consists of two neural networks, referred to as $\mathcal{N}{\mathcal{N}}_v$ and $\mathcal{N}{\mathcal{N}}_{\psi }$. The primary network, $\mathcal{N}{\mathcal{N}}_v$, is a GRU-based RNN with learnable parameters θ_v, designed to learn the history state vector v from the strain history. The auxiliary network, $\mathcal{N}{\mathcal{N}}_{\psi }$, is a fully connected neural network with learnable parameters θ_ψ, which is used to compute the Helmholtz free energy of the material. The primary network and the auxiliary network are illustrated in Figure 1B,C, respectively, and the detailed computational procedures are provided in the Appendix. In this framework, the GRU-based primary network encodes the loading history together with time increments, enabling the characterization of path-dependent and rate-dependent responses. Meanwhile, the auxiliary free-energy network, combined with physics-informed loss terms, enhances thermodynamic consistency and facilitates subsequent implementation in ABAQUS for structural simulations.

Figure 1 Schematic architecture of the physics-informed recurrent neural network surrogate model for the JC constitutive model: (A) overall architecture; (B) primary network; (C) auxiliary network

The computational procedure of the PIRNN is described as follows:

(1) Predicting the history state vector at the (n-1)-th and n-th increments with the primary network:$\begin{array}{l}{v}^{n-1},v^n=\mathcal{N}{\mathcal{N}}_v\left({\mathit{x}}^{n-L+1},\cdots ,{\mathit{x}}^k,\cdots ,{\mathit{x}}^n;{\theta }_v\right),\\ k\in \left[n-L+1,n\right],\end{array}$(9)where L is the input sequence length, x^k=(ε^k−1, Δε^k, σ^k−1, Δt^k), Δε and Δt are strain increments and time increment, respectively. In this study, the plane stress assumption is considered:$\begin{array}{l}{\mathit{\epsilon }}^{(·)}={\left[{\epsilon }_{11},{\epsilon }_{22},{\gamma }_{12}\right]}^{\text{T}(·)},\\ \mathrm{\Delta }{\mathit{\epsilon }}^{(·)}={\left[\mathrm{\Delta }{\epsilon }_{11},\mathrm{\Delta }{\epsilon }_{22},\mathrm{\Delta }{\gamma }_{12}\right]}^{\text{T}(·)},\\ {\mathit{\sigma }}^{(·)}={\left[{\sigma }_{11},{\sigma }_{22},{\sigma }_{12}\right]}^{\text{T}(·)},\end{array}$(10)where γ₁₂ = 2ε₁₂ and Δγ₁₂ = 2Δε₁₂ denote the engineering shear strain and its increment, respectively.

(2) Computing the strain and the history state vector increment at the n-th increment:${\mathit{\epsilon }}^n={\mathit{\epsilon }}^{n-1}+\mathrm{\Delta }{\mathit{\epsilon }}^n,$(11)$\mathrm{\Delta }{v}^n=v^n-v^{n-1}.$(12)

(3) Predicting the Helmholtz free energy with the auxiliary network:${\psi }^n=\mathcal{N}{\mathcal{N}}_{\psi }\left({\mathit{\epsilon }}^n,v^n;{\theta }_{\mathrm{\psi }}\right).$(13)

(4) Computing the stress and the dissipated energy via automatic differentiation:${\mathit{\sigma }}^n=\frac{\partial {\psi }^n}{\partial {\mathit{\epsilon }}^n},$(14)${\psi }_d^n=D^n\mathrm{\Delta }{t}^n\approx -{\displaystyle \sum _{i=1}^N\frac{\partial {\psi }^n}{\partial v_i^n}}\cdot \mathrm{\Delta }{v}_i^n.$(15)

During the training process of the PIRNN, minimizing the error between the predicted and exact stresses and satisfying the physical constraints specified in Equations (7, 8) are performed simultaneously. In this spirit, the loss function is defined as a weighted sum of these two components:$L_{\text{total}}={\lambda }_{\sigma }{L}_{\sigma }+{\lambda }_{\text{p}}{L}_{\text{p}},$(16)where L_σ and L_p are stress loss term and physical loss term, respectively. λ_σ and λ_p are corresponding weights of them. Each component in Eq. (16) is thus expressed by:$L_{\sigma }={\Vert \mathit{\sigma }-\mathit{\sigma }*\Vert }_2^2,$(17)$L_{\text{p}}=\text{max}\left(-\psi ,0\right)+\text{max}\left(-{\psi }_{\text{d}},0\right),$(18)where σ* is exact stress generated by the JC constitutive model. λ_σ and λ_p are updated at the (m+1)-th epoch as follows:$\begin{array}{l}{\lambda }_{\mathrm{\sigma }}^{\left(m+1\right)}=\frac{L_{\mathrm{\sigma }}^{\left(m\right)}}{L^{\left(m\right)}+L_{\mathrm{\delta }}},\\ {\lambda }_{\text{p}}^{\left(m+1\right)}=\frac{L_{\text{p}}^{\left(m\right)}}{L^{\left(m\right)}+L_{\mathrm{\delta }}},\end{array}$(19)where L_δ is a positive value to avoid the denominator being zero, and it is set to 10⁻⁵ in this study.

2.3 Integration with ABAQUS

As illustrated in Figure 2, the trained PIRNN model is incorporated as a surrogate constitutive material model into the ABAQUS solver using UMAT, cf. Lourenço et al. [32]. During structural simulations, ABAQUS acts as the macroscopic solver, while the response at each Gauss integration point is evaluated by invoking the trained PIRNN through the UMAT. The state-dependent variables (SDVs) are initialized as a zero vector, with a length consistent with the total number of inputs to the PIRNN model:$L_{\text{SDVs}}=L\times d_{\text{in}},$(20)where L_SDVs is the length of the SDVs, L denotes the length of the input sequence, and d_in is the input dimension. For each increment, the SDVs, together with the current strain increment Δεⁿ and time increment Δtⁿ, are assembled as inputs to the PIRNN to predict the stress, after which the updated SDVs and stress are returned to ABAQUS. This integration with ABAQUS is used to evaluate the performance of the PIRNN for structural simulations.

Figure 2 Flowchart of the PIRNN integration with ABAQUS

3.1 Dataset generation

Synthetic stress-strain data are first generated on a single material point from multiaxial loading with the conventional JC constitutive model implemented in ABAQUS. The constitutive parameters of the JC model employed in this study are listed in Table 1, which were calibrated by Quan et al. [45] based on experimental data reported by Shi et al. [46].

The temperature is fixed at 293 K, and 10,000 random-walk strain paths are prescribed. Each strain path consists of ten multiaxial loading or reverse loading steps with durations ranging from 5×10⁻⁸ to 8×10⁻⁷ s. Each step is further discretized into 20 equal increments. For each time increment, the strain and stress components, together with the time increment, are extracted to construct the training datasets. Figure 3 illustrates the randomly selected 50 strain paths along with their corresponding stress responses, with one path highlighted in red. These strains, strain increments, and stresses are normalized to the range of −1 to 1, which facilitates the training of the PIRNN.

Figure 3 Synthetic stress-strain data for PIRNN training: (A) strain paths; (B) corresponding stress responses

3.2 Validations

Table 2 summarizes the hyperparameters for the PIRNN model. Stress-strain datasets are divided into three subsets with a ratio of 3:1:1 for training, validation, and testing of the PIRNN model, respectively. Figure 4 presents the training loss, validation loss, and adaptive weights of the PIRNN across training epochs. Within the first 15 epochs, both losses decrease steadily. Subsequently, the validation loss begins to fluctuate and shows an increasing trend after 30 epochs, indicating the onset of overfitting. Therefore, the training is terminated, and the checkpoint corresponding to the minimum validation loss is selected as the optimal weights of the PIRNN.

Figure 4 Learning curves for the PIRNN model: (A) training loss; (B) validation loss; (C) adaptive weights

The stress components directly predicted by the trained PIRNN are compared with the exact stress components on the unseen test dataset, as shown in Figure 5. The mean absolute errors (MAE) of the stress components σ₁₁, σ₂₂, and σ₁₂ are 4.57, 5.03, and 2.33 MPa, respectively, while the corresponding coefficients of determination R² reach 0.9999. All scatter points are distributed closely around the diagonal line. These results demonstrate that the PIRNN achieves high predictive accuracy on the test datasets.

Figure 5 Comparison of the exact and predicted stress components by PIRNN on the test dataset: (A) σ₁₁; (B) σ₂₂; (C) σ₁₂

The trained PIRNN model is integrated into ABAQUS following the procedure illustrated in Figure 2. A single material point is then considered along a strain path selected from the test datasets, as shown in Figure 6A−C. The rates of the three strain components ε₁₁, ε₂₂, and γ₁₂ vary at each step. Figure 6D−F compare stress components obtained from the PIRNN and JC constitutive models. The R² of σ₁₁, σ₂₂, and σ₁₂ are 0.9981, 0.9987, and 0.9959, while the MAE of σ₁₁, σ₂₂, and σ₁₂ are 24.12, 9.24, and 10.02 MPa. This indicates that the PIRNN has successfully captured the rate-dependent and path-dependent behavior of the JC constitutive model.

Figure 6 Prescribed strain paths and comparison of stress-time responses between the PIRNN and JC constitutive model: (A–C) prescribed strain paths of ε₁₁, ε₂₂, and γ₁₂; (D–F) comparison of the predicted and reference stress responses of σ₁₁, σ₂₂, and σ₁₂

3.3 Experimental validations

A single-element model is constructed in ABAQUS, where all stress components except σ₁₁ are constrained to 0 to simulate uniaxial tensile loading under four different strain rates, namely, 3×10³, 5×10³, 8×10³ and 1.1×10⁴ s⁻¹. The PIRNN-generated results are compared against the experimental data reported by Shi et al. [46], as shown in Figure 7. We note that the PIRNN predictions are in excellent agreement with the experimental values under all strain-rate conditions.

Figure 7 Comparison of PIRNN-predicted stress-strain responses against experimental data under tensile loading with different loading rates: (A) $\dot{\epsilon }$ = 3×10³ s⁻¹ and $\dot{\epsilon }$ = 5×10³ s⁻¹; (B) $\dot{\epsilon }$ = 8×10³ s⁻¹ and $\dot{\epsilon }$ = 1.1×10⁴ s⁻¹