A Physics-Guided Machine Learning Model for Predicting Viscoelasticity of Solids at Large Deformation

Abstract

Physics-guided machine learning (PGML) methods are emerging as valuable tools for modelling the constitutive relations of solids due to their ability to integrate both data and physical knowledge. While various PGML approaches have successfully modeled time-independent elasticity and plasticity, viscoelasticity remains less addressed due to its dependence on both time and loading paths. Moreover, many existing methods require large datasets from experiments or physics-based simulations to effectively predict constitutive relations, and they may struggle to model viscoelasticity accurately when experimental data are scarce. This paper aims to develop a physics-guided recurrent neural network (RNN) model to predict the viscoelastic behavior of solids at large deformations with limited experimental data. The proposed model, based on a combination of gated recurrent units (GRU) and feedforward neural networks (FNN), utilizes both time and stretch (or strain) sequences as inputs, allowing it to predict stress dependent on time and loading paths. Additionally, the paper introduces a physics-guided initialization approach for GRU–FNN parameters, using numerical stress–stretch data from the generalized Maxwell model for viscoelastic VHB polymers. This initialization is performed prior to training with experimental data, helping to overcome challenges associated with data scarcity.

1. Introduction

Machine learning has made significant strides in various fields of science and engineering, driven by the proliferation of digital data, increasing computing power, and advanced algorithms. It has achieved notable success in applications such as speech recognition [1], image classification [2], and cognitive science [3]. Recently, machine learning has emerged as a promising tool for predicting deformation patterns [4,5] by modeling the mechanical constitutive behavior of solids, particularly the relationship between mechanical stress and strain. However, modeling both time and history-dependent nonlinear constitutive relations such as viscoelasticity presents challenges for two reasons. First, many state-of-the-art machine learning techniques (e.g., feedforward/convolutional/recurrent neural networks, FNN/CNN/RNN) [6,7,8] require large datasets and often lack robustness, failing to guarantee convergence when only limited experimental stress–strain data are available from constitutive experiments, such as uniaxial tensile or bending tests. Second, modelling time and history-dependent nonlinear stress–strain relations can be difficult for many machine learning techniques, including the popular FNN, due to their limited capacity to handle sequential information. Therefore, advancing machine learning approaches to effectively model viscoelasticity in solids at large deformation is crucial for furthering the study of solid mechanics problems.

Physics-guided machine learning (PGML) methods, which integrate established physical principles with data-driven approaches, have proven beneficial in solid mechanics by simultaneously incorporating data information and physical knowledge [9,10,11]. These studies can be broadly classified into three categories: (i) discovering constitutive equations from data [12,13,14,15,16,17]; (ii) solving governing equations [18,19,20,21,22]; (iii) predicting data-driven constitutive relations [15,23,24,25,26,27,28]. Here, we focus on elaborating the third category. To predict data-driven constitutive relations, significant efforts have been made using deep neural networks, including FNN, temporal convolutional network (TCN), and RNN, among others. For example, Linka et al. [15] proposed a data-driven constitutive neural network model for predicting the mechanical constitutive behavior of hyperelastic materials. This model integrates physical laws and the symmetry of mechanical properties into the network structure, enabling the simultaneous use of information from three sources: stress–strain data, theoretical knowledge from materials science, and additional data such as microstructural or processing information; Danoun et al. [23] developed a thermodynamically consistent RNN model to simulate the mechanical response of elastoplastic materials under multi-axial and non-proportional loading conditions. Their work demonstrates that incorporating thermodynamic consistency can significantly enhance the predictive capabilities of such surrogate models; Wang et al. [24] developed a TCN model for materials with an ultra-long-history-dependent stress–strain relation. Their TCN model outperforms RNN models, including Long Short-Term Memory (LSTM) and Gated Recurrent Unit (GRU) models, in terms of both performance and speed for various sequential tasks. Abdolazizi et al. [28] proposed viscoelastic Constitutive Artificial Neural Networks (vCANNs), a novel physics-informed machine learning framework for anisotropic nonlinear viscoelasticity at finite strains. This approach is based on generalized Maxwell models, enhanced with neural networks to capture nonlinear strain- and strain rate-dependent properties. However, most of these models require a moderate to large amount of labeled data from experiments or physics-based simulations to accurately learn time-independent constitutive relations, such as those for path-independent elastic materials and path-dependent plastic materials [26,27]. They often face challenge in modelling the viscoelasticity of solids, which depends on both time and loading paths, especially when experimental data are scarce. Additionally, incorporating physical knowledge of viscoelasticity, such as viscous dissipation energy, directly into the training process is challenging. This difficulty arises because instantaneous viscous deformation and its conjugated driving force, two key components for calculating viscous dissipation energy, are difficult to obtain through experiments.

Therefore, this paper aims to develop a physics-guided machine learning model for predicting the viscoelasticity of solids at large deformation with scarce experimental data. First, we propose a surrogate model of viscoelasticity based on GRU and FNN (GRU–FNN). This model inputs both time and stretch (or strain) sequences, enabling it to predict stress based on time and loading path. Second, we implement a physics-guided initialization for the GRU–FNN parameters by training on numerical stress–stretch data derived from the generalized Maxwell model for viscoelastic solids under large deformations. This approach allows the model to predict the viscoelastic behavior of solids even with scarce experimental stress–stretch data. Unlike existing data-driven models that use residuals of the governing equations or variational forms as penalizing terms to restrict the solution space, our strategy effectively avoids the need for direct measurements of viscous deformation.

The remainder of this paper is organized as follows. In Section 2, a continuum theoretical framework for the viscoelasticity of solids at large deformation is given and constitutive equations including the state and evolving equations are derived. In Section 3, a GRU–FNN model for viscoelasticity is established. In Section 4, uniaxial viscoelasticity of commercially available dielectric polymers VHB4905 is modelled. Finally, conclusions are given in Section 5.

2. Thermodynamic Formulation of Constitutive Laws for Viscoelastic Solids at Large Deformation

Consider a solid $B_0$ bounded by the surface $\mathit{\partial }{B}_0$, defined in a fixed reference configuration, which deforms into $B_t$ with the surface $\mathit{\partial }{B}_t$. The deformation gradient is then given by

$$\mathit{F}={\nabla }_R\mathit{x}$$

(1)

Here, $\mathit{x}=\mathit{\chi }\left(\mathit{X},t\right)$ is a function that maps an arbitrary material point $\mathit{X}$ inside $B_0$ into a spatial point $\mathit{x}$ inside $B_t$, and ${\nabla }_R$ represents the gradient operator with respect to the coordinates $\mathit{X}$. To describe the viscoelasticity of solids under large deformation, a generalized Maxwell model [29,30,31,32,33] is employed, as shown in Figure 1. In this model, the solid is represented by an equilibrium spring, characterized by the deformation gradient $\mathit{F}$, along with $n$ parallel Maxwell elements (each consisting of a spring and a dashpot in series), characterized by the elastic deformation gradient ${\mathit{F}}_{\mathit{i}}^{\mathit{e}}$ and the viscous deformation gradient ${\mathit{F}}_{\mathit{i}}^{\mathit{v}}$, leading to the following relation:

$$\mathit{F}={\mathit{F}}_i^e\cdot {\mathit{F}}_i^v, i=1,\dots ,n$$

(2)

Accordingly, the right Cauchy–Green deformation tensors, $\mathit{C}={\mathit{F}}^T\cdot \mathit{F}$ and ${\mathit{C}}_i^e={\left({\mathit{F}}_i^e\right)}^T\cdot {\mathit{F}}_i^e$, where the superscript $T$ denotes the transpose, are used to measure the deformation of the equilibrium spring and the elastic deformation of the Maxwell elements, respectively.

Neglecting inertial effects, the balance laws of force and moment in the current configuration are given by

$$\mathit{\sigma }\cdot \nabla +\mathit{b}=\mathbf{0}, \mathit{\sigma }={\mathit{\sigma }}^T$$

(3)

where $\mathit{\sigma }$ is the Cauchy stress tensor and $\mathit{b}$ is the body force per unit volume in $B_t$, which can be expressed in the reference configuration as

$$\mathit{P}\cdot {\nabla }_R+{\mathit{b}}_R=\mathbf{0}, P\cdot F^T=\mathit{F}\cdot {\mathit{P}}^T$$

(4)

where $\mathit{P}$ denotes the first P–K stress tensor and ${\mathit{b}}_R$ denotes the body force per unit volume in $B_0$. Let $j$ denote the determinant of $\mathit{F}$, and we have the relations $\mathit{P}=j\mathit{\sigma }\cdot {\left({\mathit{F}}^T\right)}^{-1}$ and ${\mathit{b}}_R=j\mathit{b}$.

Let $\epsilon$ and $q$ represent the internal energy density and the heat source per unit volume in $B_0$, respectively, and ${\mathit{j}}^q$ denote the heat flux per unit area on $\mathit{\partial }{B}_t$. The energy balance law in the current configuration is then expressed as

$$\dot{\epsilon }=\mathit{\sigma }:\nabla \dot{\mathit{\chi }}-\nabla ·{\mathit{j}}^q+q$$

(5)

whose corresponding form in the reference configuration is

$${\dot{\epsilon }}_R=\mathit{P}:\dot{\mathit{F}}-{\nabla }_R\cdot {\mathit{j}}_R^q+q_R$$

(6)

where ${\displaystyle \frac{d}{dt}}$ or equivalently a superposed dot represents the material time derivative. ${\epsilon }_R$ and $q_R$ denote the internal energy density and the heat source per unit volume in $B_0$, respectively, while ${\mathit{j}}_R^q$ represents the heat flux per unit area on $\mathit{\partial }{B}_0$. The relationships ${\epsilon }_R=j\epsilon$, $q_R=jq$ and ${\mathit{j}}_R^q=j{\mathit{j}}^q\cdot {\left({\mathit{F}}^T\right)}^{-1}$ hold for the above two equations. The first term on the right-hand side represents the mechanical work, while the second and third terms account for the energy from heat flow across the surface and the heat source inside the body, respectively.

Let $\eta$ denote the entropy per unit volume in $B_t$ and $\vartheta$ the absolute temperature. The entropy inequality in the current configuration is expressed locally as

$$\dot{\eta }\ge -\nabla \cdot {\displaystyle \frac{{\mathit{j}}^q}{\vartheta }}+{\displaystyle \frac{q}{\vartheta }}$$

(7)

which can also be written in the reference configuration as

$${\dot{\eta }}_R\ge -{\nabla }_R\cdot {\displaystyle \frac{{\mathit{j}}_R^q}{\vartheta }}+{\displaystyle \frac{q_R}{\vartheta }}$$

(8)

where ${\eta }_R$ is the entropy per unit volume in $B_0$, related to $\eta$ by ${\eta }_R=j\eta$.

Introducing the Helmholtz free energy density $\phi =\epsilon -\vartheta \eta$ and considering the energy balance (5), the inequality (7) becomes

$$\mathit{\sigma }:\nabla \dot{\mathit{\chi }}-\eta \dot{\vartheta }-\dot{\phi }-{\displaystyle \frac{1}{\vartheta }}{\mathit{j}}^q\cdot \nabla \vartheta \ge 0$$

(9)

Similarly, by introducing the Helmholtz free energy density ${\phi }_R={\epsilon }_R-\vartheta {\eta }_R$, where ${\phi }_R=j\phi$, and considering the energy balance (6), the inequality becomes

$$\mathit{P}:\dot{\mathit{F}}-{\eta }_R\dot{\vartheta }-{\dot{\phi }}_R-{\displaystyle \frac{1}{\vartheta }}{\mathit{j}}_R^q\cdot {\nabla }_R\vartheta \ge 0$$

(10)

This imposes the thermodynamic constraint on solids. For convenience, we will use the formulations provided in the reference configuration in the following sections.

Considering the thermo-viscoelastic effects in solids, the Helmholtz free energy density ${\phi }_R$ can be assumed to be a function of variables $\left(\mathit{C},{\mathit{C}}_i^e,\vartheta \right)$, i.e.,

$${\phi }_R={\phi }_R\left(\mathit{C},{\mathit{C}}_i^e,\vartheta \right)$$

(11)

whose material time derivative can be further written as

$${\dot{\phi }}_R={\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }\mathit{C}}}:\dot{\mathit{C}}+{\displaystyle \sum _{i=1}^n}{\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}:{\dot{\mathit{C}}}_i^e+{\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }\vartheta }}\dot{\vartheta }$$

(12)

The first and second terms on the right-hand side of Equation (12) can be rewritten as

$$\begin{array}{l}{\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }\mathit{C}}}:\dot{\mathit{C}}=2\mathit{F}\cdot {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }\mathit{C}}}:\dot{\mathit{F}}\\ {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}:{\dot{\mathit{C}}}_i^e={\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}:{\displaystyle \frac{\mathit{\partial }{\mathit{C}}_i^e}{\mathit{\partial }\mathit{F}}}:\dot{\mathit{F}}+{\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}:{\displaystyle \frac{\mathit{\partial }{\mathit{C}}_i^e}{\mathit{\partial }{\mathit{F}}_i^v}}:{\dot{\mathit{F}}}_i^v=2{\mathit{F}}_i^e\cdot {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}\cdot {\left({\mathit{F}}_i^{vT}\right)}^{-1}:\dot{\mathit{F}}-2{\mathit{C}}_i^e\cdot {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}:{\mathit{D}}_i^v\end{array}$$

(13)

with

$${\mathit{D}}_i^v={\displaystyle \frac{1}{2}}\left({\mathit{L}}_i^v+{\mathit{L}}_i^{vT}\right), L_i^v={\dot{\mathit{F}}}_i^v\cdot {\left({\mathit{F}}_i^v\right)}^{-1}$$

(14)

where the superscript ‘$-1$’ denotes the inverse of a tensor. Then, substitution of Equations (12) and (13) into Equation (10) yields

$$\left(\mathit{P}-2\mathit{F}\cdot {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }\mathit{C}}}-2{\mathit{F}}_i^e\cdot {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}\cdot {\left({\mathit{F}}_i^{vT}\right)}^{-1}\right):\dot{\mathit{F}}-\left({\eta }_R+{\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }\vartheta }}\right)\dot{\vartheta }-{\displaystyle \frac{1}{\vartheta }}{j}_R^q\cdot {\nabla }_R\vartheta +{{\displaystyle \sum }}_{i=1}^{n}2{\mathit{C}}_i^e\cdot {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}:{\mathit{D}}_i^v\ge 0$$

(15)

In the case where $\mathit{P}$ and ${\eta }_R$ are independent of $\dot{\mathit{F}}$ and $\dot{\vartheta }$, the first two terms of the above inequality must vanish, leading to the following constitutive relations:

$$P=2\mathit{F}\cdot {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }\mathit{C}}}+2{\mathit{F}}_i^e\cdot {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}\cdot {\left({\mathit{F}}_i^{vT}\right)}^{-1}$$

(16)

$${\eta }_R=-{\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }\vartheta }}$$

(17)

Thus, the inequality (15) reduces to

$$-{\displaystyle \frac{1}{\vartheta }}{\mathit{j}}_R^q\cdot {\nabla }_R\vartheta +{\displaystyle \sum _{i=1}^n}{\mathit{M}}_i^{neq}:{\mathit{D}}_i^v\ge 0$$

(18)

where ${\mathit{M}}_i^{neq}$ is the non-equilibrium Mandel stress tensor [34], defined as

$${\mathit{M}}_i^{neq}=2{\mathit{C}}_i^e\cdot {\displaystyle \frac{\mathit{\partial }{\phi }_R}{\mathit{\partial }{\mathit{C}}_i^e}}$$

(19)

More specially, to satisfy the thermodynamic constraint imposed by inequality (18), the following constitutive equations are derived:

$${\mathit{j}}_R^q=-\mathit{Y}\cdot {\nabla }_R\vartheta$$

(20)

$${\mathit{D}}_i^v={\mathit{Q}}_i:{\mathit{M}}_i^{neq}$$

(21)

where $\mathit{Y}$ is a second order tensor and ${\mathit{Q}}_i$ is a fourth order tensor, both of which are positive-definite. Here, Equation (20) represents the Fourier heat conduction law, and Equation (21) describes the rheological viscous flow rule [30,34].

The Helmholtz free energy density is assumed to be additively decomposed as follows:

$${\phi }_R={\phi }_R^{eq}+{\phi }_R^{neq}$$

(22)

where ${\phi }_R^{eq}$ and ${\phi }_R^{neq}$ represent the equilibrium and non-equilibrium Helmholtz free energy densities, corresponding to the stretching of the single spring and the springs in Maxwell element, respectively. We adopt the Gent model [35] for both ${\phi }_R^{eq}$ and ${\phi }_R^{neq}$ to account for the strain-stiffening effect, where solids may stiffen sharply as the stretch approaches their extension limit [36], as follows:

$${\phi }_R^{eq}=-{\displaystyle \frac{G^{eq}L}{2}}\mathit{ln}\left({\displaystyle \frac{L-tr\mathit{C}+3}{L}}\right)$$

(23)

$${\phi }_R^{neq}=-\sum _{i=1}^n{\displaystyle \frac{G_i^{neq}{L}_i}{2}}\mathit{ln}\left({\displaystyle \frac{L_i-tr{\mathit{C}}_i^e+3}{L_i}}\right)$$

(24)

where the symbol ‘$tr$’ denotes the trace of a tensor, $G^{eq}$ and $G_i^{neq}$ represent the equilibrium modulus of the single spring and the nonequilibrium modulus of the ith Maxwell element, respectively, $L$ and $L_i$ denote the extension limits of the single spring and the ith Maxwell element spring, respectively. Substituting Equations (22) into (16) and (19), we have

$$\mathit{P}={\displaystyle \frac{G^{eq}L}{L-tr\mathit{C}+3}}\mathit{F}+\sum _{i=1}^n{\displaystyle \frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}}{\mathit{F}}_i^e\cdot {({\mathit{F}}_i^{vT})}^{-1}$$

(25)

$${\mathit{M}}_i^{neq}={\displaystyle \frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}}{\mathit{C}}_i^e$$

(26)

Next, consider a case that the elastic deformations of springs are incompressible, which generally applies to polymers [32]. Lagrange multipliers $\mathsf{\Pi }$ and ${\mathsf{\Pi }}_i$ are introduced to enforce these constraint conditions, modifying the free energy density as

$${\tilde{\phi }}_R={\phi }_R-\mathsf{\Pi }\left(j-1\right)-\sum _{i=1}^n{\mathsf{\Pi }}_i\left(j_i^e-1\right)$$

(27)

where $j_i^e$ is the determinant of the elastic deformation gradient ${\mathit{F}}_i^e$. Replacing ${\phi }_R$ with ${\tilde{\phi }}_R$ in Equations (16) and (19), we can rewrite $\mathit{P}$ and ${\mathit{M}}_i^{neq}$ as

$$\mathit{P}={\displaystyle \frac{G^{eq}L}{L-tr\mathit{C}+3}}\mathit{F}+{\sum }_{i=1}^n{\displaystyle \frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}}{\mathit{F}}_i^e\cdot {({\mathit{F}}_i^{vT})}^{-1}-j\mathsf{\Pi }{\left({\mathit{F}}^T\right)}^{-1}-\sum _{i=1}^n{j}_i^e{\mathsf{\Pi }}_i{\left({\mathit{F}}^T\right)}^{-1}$$

(28)

$${\mathit{M}}_i^{neq}={\displaystyle \frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}}{\mathit{C}}_i^e-j_i^e{\mathsf{\Pi }}_i\mathit{I}$$

(29)

where $\mathit{I}$ is the second-order unit tensor. Furthermore, the Cauchy stress tensor $\mathit{\sigma }$ can be obtained using the relations $\mathbf{\sigma }={\displaystyle \frac{1}{j}}\mathit{P}\cdot {\mathit{F}}^T$, and $j=j_i^e=1$ in Equation (28), as follows:

$$\mathit{\sigma }={\displaystyle \frac{G^{eq}L}{L-tr\mathit{C}+3}}\mathit{F}\cdot {\mathit{F}}^T+{\sum }_{i=1}^n{\displaystyle \frac{G_i^{neq}{L}_i}{L_i-tr{\mathit{C}}_i^e+3}}{\mathit{F}}_i^e\cdot {\mathit{F}}_i^{eT}-\left(\mathsf{\Pi }+\sum _{i=1}^n{\mathsf{\Pi }}_i\right)\mathit{I}$$

(30)

By employing $j=j_i^e=1$ and the multiplicative decomposition of $\mathit{F}$ from Equation (2), the condition of viscous incompressibility, $tr{\mathit{D}}_i^v=\mathbf{0}$ or $j_i^v=1$, where $j_i^v$ is the determinant of ${\mathit{F}}_i^v$, can be deduced. From this, the fourth-order tensor $Q_i$ in Equation (21) can be expressed as [32]

$${\mathit{Q}}_i={\displaystyle \frac{1}{2v_i}}\left({\mathit{I}}_4-{\displaystyle \frac{1}{3}}\mathit{I}\otimes \mathit{I}\right)$$

(31)

where ${\mathit{I}}_4$ is the fourth-order unit tensor and $v_i$ is the viscosity of the $i\mathrm{th}$ subnetwork. The relaxation time for the $i\mathrm{th}$ Maxwell element is then defined as ${\tau }_i={\displaystyle \frac{v_i}{G_i^{neq}}}$ [31].

The corresponding deformation gradient in tensile tests, without shear deformation, can be expressed as

$$\begin{array}{l}\mathit{F}=\left[\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right], {\mathbf{F}}_i^e=\left[\begin{array}{ccc}{\lambda }_{i1}^e& 0& 0\\ 0& {\lambda }_{i2}^e& 0\\ 0& 0& {\lambda }_{i3}^e\end{array}\right], {\mathbf{F}}_i^v=\left[\begin{array}{ccc}{\lambda }_{i1}^v& 0& 0\\ 0& {\lambda }_{i2}^v& 0\\ 0& 0& {\lambda }_{i3}^v\end{array}\right]\\ {\mathit{D}}_i^v=\left[\begin{array}{ccc}{\dot{\lambda }}_{i1}^v/{\lambda }_{i1}^v& 0& 0\\ 0& {\dot{\lambda }}_{i2}^v/{\lambda }_{i2}^v& 0\\ 0& 0& {\dot{\lambda }}_{i3}^v/{\lambda }_{i3}^v\end{array}\right]\end{array}$$

(32)

Here ${\lambda }_1$, ${\lambda }_2$, ${\lambda }_3$ are the principal stretches of the deformation gradient $\mathit{F}$; ${\lambda }_{i1}^e$, ${\lambda }_{i2}^e$, ${\lambda }_{i3}^e$ are the principle stretches of the elastic deformation gradient ${\mathit{F}}_i^e$; and ${\lambda }_{i1}^v$, ${\lambda }_{i2}^v$, ${\lambda }_{i3}^v$ are the principal stretches of the viscous deformation gradient ${\mathit{F}}_i^v$. Additionally, considering the equal lateral stretches during the uniaxial tensile test and the incompressible deformation condition $j=j_i^e=j_i^v=1$, we have

$${\lambda }_1=\lambda , {\lambda }_2={\lambda }_3={\left(\lambda \right)}^{-{\displaystyle \frac{1}{2}}}$$

$${\lambda }_1={\lambda }_{i1}^e{\lambda }_{i1}^v=\lambda , {\lambda }_{i2}^e{\lambda }_{i2}^v={\lambda }_{i3}^e{\lambda }_{i3}^v={\left(\lambda \right)}^{-{\displaystyle \frac{1}{2}}}$$

$${\lambda }_{i1}^e={\lambda }_i^e, {\lambda }_{i2}^e={\lambda }_{i3}^e={\left({\lambda }_i^e\right)}^{-{\displaystyle \frac{1}{2}}}, {\lambda }_{i1}^v={\lambda }_i^v, {\lambda }_{i2}^v={\lambda }_{i3}^v={\left({\lambda }_i^v\right)}^{-{\displaystyle \frac{1}{2}}}$$

(33)

with $\lambda$, ${\lambda }_i^e$, ${\lambda }_i^v$ being the principal stretches of the deformation gradients $\mathit{F}$, ${\mathit{F}}_i^e$, ${\mathit{F}}_i^v$ along the tensile direction, respectively.

Let $P_1$, $P_2$ $P_3$ denote the components of $P$ along the three principal directions, respectively. Substituting Equations (32) and (33) into Equation (28), we obtain

$$\begin{array}{l}{P}_1={\displaystyle \frac{G^{eq}L\lambda }{L-2{\lambda }^{-1}-{\lambda }^2+3}}+{\sum }_{i=1}^n{\displaystyle \frac{G_i^{neq}{L}_i\lambda {\left({\lambda }_i^v\right)}^{-2}}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}}-{\displaystyle \frac{\mathsf{\Pi }+{\sum }_{i=1}^n{\mathsf{\Pi }}_i}{\lambda }}\\ P_2=P_3={\displaystyle \frac{G^{eq}L{\lambda }^{-\frac{1}{2}}}{L-2{\lambda }^{-1}-{\lambda }^2+3}}+{\sum }_{i=1}^n{\displaystyle \frac{G_i^{neq}{L}_i{\lambda }^{-\frac{1}{2}}{\lambda }_i^v}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}}-\left(\mathsf{\Pi }+{\displaystyle \sum _{i=1}^n}{\mathsf{\Pi }}_i\right){\lambda }^{{\displaystyle \frac{1}{2}}}\end{array}$$

(34)

According to the boundary condition $P_2=P_3=0$ for uniaxial loading–unloading tests, we can further derive

$$\mathsf{\Pi }+\sum _{i=1}^n{\mathsf{\Pi }}_i={\displaystyle \frac{G^{eq}L{\lambda }^{-1}}{L-2{\lambda }^{-1}-{\lambda }^2+3}}+{\sum }_{i=1}^n{\displaystyle \frac{G_i^{neq}{L}_i{\lambda }^{-1}{\lambda }_i^v}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}}$$

(35)

and

$$P_1={\displaystyle \frac{G^{eq}L\left(\lambda -{\lambda }^{-2}\right)}{L-2{\lambda }^{-1}-{\lambda }^2+3}}+{\sum }_{i=1}^n{\displaystyle \frac{G_i^{neq}{L}_i\left[\lambda {\left({\lambda }_i^v\right)}^{-2}-{\lambda }^{-2}{\lambda }_i^v\right]}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}}$$

(36)

Here, the effect of deformation incompressibility on the stress is considered, and the Lagrange multipliers are eliminated according to the boundary conditions.

Let $M_{i1}^{neq}$, $M_{i2}^{neq}$ $M_{i3}^{neq}$ denote the components of ${\mathit{M}}_i^{neq}$ along the three principal directions, respectively. Similarly, substituting Equations (32) and (33) into Equation (29), we obtain

$$\begin{array}{l}{M}_{i1}^{neq}={\displaystyle \frac{G_i^{neq}{L}_i{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}}-{\mathsf{\Pi }}_i\\ M_{i2}^{neq}=M_{i3}^{neq}={\displaystyle \frac{G_i^{neq}{L}_i{\lambda }^{-1}{\lambda }_i^v}{L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3}}-{\mathsf{\Pi }}_i\end{array}$$

(37)

Then, substitution of Equations (37), (32) and (31) into Equation (21) yields

$${\displaystyle \frac{{\dot{\lambda }}_i^v}{{\lambda }_i^v}}={\displaystyle \frac{L_i}{3{\tau }_i\left[L_i-2{\lambda }^{-1}{\lambda }_i^v-{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}+3\right]}}\left[{\lambda }^2{\left({\lambda }_i^v\right)}^{-2}-{\lambda }^{-1}{\lambda }_i^v\right]$$

(38)

This describes the viscous flow in solids subjected to uniaxial tension. It is important to note that the elastic incompressibility of the springs does not affect the viscous flow, as only the deviatoric stress, excluding the hydrostatic pressure ${\mathsf{\Pi }}_i$, derives the viscous flow.

Under the initial conditions ${\left.P_1\right|}_{t=0}=0$, ${\left.\lambda \right|}_{t=0}={\left.{\lambda }_i^v\right|}_{t=0}=1$, the coupled Equations (36) and (38) can be solved simultaneously to predict the mechanical behaviors based on the given material parameters. A finite difference approach is used to obtain the numerical results for the evolving stress $P_1$ as a function of the stretch $\lambda$ at ten different stretching rates. The model parameters are obtained by fitting Equations (36) and (38) with the experimental data at 296K as n = 3, $L=L_i=155 \left(i=1, 2, 3\right)$, $G^{eq}=15.12 \mathrm{kPa}$, $G_1^{neq}=16.83 \mathrm{kPa}$, $G_2^{neq}=27.39 \mathrm{kPa}$, $G_3^{neq}=41.04 \mathrm{kPa}$, ${\tau }_1=413.32 \mathrm{s}$, ${\tau }_2=5.43\mathrm{s}$, ${\tau }_3=1.65 \mathrm{s}$. Great consistency between the experimental data and the theoretical model is shown in Figure 2. With these parameters determined, the theoretical data at different stretching rates are calculated as shown in Figure 3. The stress can be observed to increase with stretch during the loading phase and decrease with stretch during the unloading phase. Each loading–unloading cycle exhibits characteristic viscoelastic behavior, with the peak stress magnitude rising as the stretching rate increases.

The initialization of parameters of the following machine learning model will be obtained by training these theoretical results, which is proven to be an effective strategy for predicting viscoelasticity of solids with scarce experimental data. This strategy avoids directly imposing physical constraints on the training procedure since the specific governing equations are usually uncertain and the physical quantities such as viscous deformation and its conjugated driving force are hard to obtain via experiments. In the next section, a machine learning model for predicting viscoelasticity will be introduced.

3. Machine Learning Method for Predicting Viscoelasticity

As mentioned, the FNN architecture is not well suited for handling history-dependent behaviors because there is no direct relationship between the network inputs at a current time step and the outputs from previous time steps. To address this issue, the RNN architecture is introduced as a reliable model for managing this history dependency. RNN training utilizes the Backpropagation Through Time (BPTT) algorithm, which applies backpropagation [37] to a time sequence. However, when dealing with long sequence inputs, RNN training may suffer from the common problems of vanishing and exploding gradients during backpropagation. This occurs because RNN builds very deep computational graphs by repeatedly applying the same operation at each time step of a long sequence. The LSTM and GRU are specialized recurrent neural network architectures designed to tackle these vanishing and exploding gradient issues. Memory cells in LSTM and GRU networks are equipped with gated units that dynamically control the flow of information, allowing the networks to “forget” old and unnecessary information and avoid the problems associated with multiplying large sequences of numbers during temporal backpropagation. GRU have the advantage of avoiding overfitting and reducing training time compared to LSTM, as GRU with two gates have significantly few parameters than LSTM with four gates.

To predict the time-dependent stress of viscoelastic materials, a surrogate model based on GRU and FNN (GRU–FNN) is employed, as illustrated in Figure 4a. At each time step, the GRU cell manages the information flow using a reset gate $r_t$, which regulates the integration of new input with the previous memory; and an update gate $Z_t$, which determines how much of the previous memory should be retained; and a hidden state $h_t$, which transfers information from the candidate hidden state ${\tilde{h}}_t$ forward, as depicted in Figure 4b. These two gates control the information flow between the long-term hidden state and the predictions at each time step. FNN then processes the final hidden state of the GRU as unidirectional flow signals, moving from the input layer through hidden layers to the output layer, as shown in Figure 4c. A detailed introduction to GRU and FNN will follow.

3.1. Gated Recurrent Unit

The reset gate output is [8,38]

$${\mathit{r}}_t=\sigma \left({\mathit{W}}^r{\mathit{x}}_t+{\mathit{U}}^r{\mathit{h}}_{t-1}+{\mathit{b}}^r\right)$$

(39)

Here, the weight matrix ${\mathit{W}}^r$ connects the input at a specific time $t$ to the current hidden layer within the reset gate, while the weight matrix ${\mathit{U}}^r$ represents the recurrent connection between the previous and current hidden layers within the reset gate. The initial hidden state is initialized the null vector ${\mathit{h}}_0=\mathbf{0}$. ${\mathit{b}}^r$ denotes the bias vector, and $\sigma$ refers to the activation function.

The update gate output is

$${\mathit{Z}}_t=\sigma \left({\mathit{W}}^z{\mathit{x}}_t+{\mathit{U}}^z{\mathit{h}}_{t-1}+{\mathit{b}}^z\right)$$

(40)

where ${\mathit{W}}^z$, ${\mathit{U}}^z$, ${\mathit{b}}^z$ are respectively two weights matrices and a bias vector regulating the update mechanism in this gate.

The candidate hidden state ${\tilde{\mathit{h}}}_t$ is

$${\tilde{\mathit{h}}}_t=\mathit{tanh}\left[{\mathit{W}}^h{\mathit{x}}_t+{\mathit{U}}^h\left({\mathit{r}}_t\ast {\mathit{h}}_{t-1}\right)+{\mathit{b}}^h\right]$$

(41)

where ${\mathit{W}}^h$ and ${\mathit{U}}^h$ are two different weight matrices, ${\mathit{b}}^h$ is a bias vector, the symbol $\ast$ denote the Hadamard product of two vectors with identical dimensions, $tanh$ is the hyperbolic tangent function.

The GRU output current hidden state ${\mathit{h}}_t$ is

$${\mathit{h}}_t=\left(\mathbf{1}-{\mathit{Z}}_t\right)\ast {\mathit{h}}_{t-1}+{\mathit{Z}}_t\ast {\tilde{\mathit{h}}}_t$$

(42)

where $\mathbf{1}$ is a unit vector, i.e., a vector filled with ones.

The total number of parameters $N_P$ in a GRU cell is related to the dimension $N_I$ of the input vector ${\mathit{x}}_t$ and the number of neurons $N_H$ in the hidden layers of the GRU cell, given by

$$N_P=3N_H\left(N_I+N_H+2\right)$$

(43)

Note that the above bias vector is considered to consist of two different bias vectors related to two different weight matrices, respectively.

3.2. Feedforward Neural Network

The FNN can be described by the following compositional function [39]:

$$O\left({\mathit{h}}_t,\mathit{W},\mathit{b}\right)={\sigma }^{\Lambda }\circ {\sigma }^{\Lambda -1}\circ \cdots \circ {\sigma }^1$$

(44)

where the symbol $\circ$ denotes the composition operator, the superscript $\Lambda$ represents the total number of the neural network, and $\sigma$ is the activation function with the following form:

$${\sigma }^{\lambda }\left({\mathit{L}}^{\lambda -1}\right)={\sigma }^{\lambda }\left({\mathit{W}}^{\lambda }{\mathit{L}}^{\lambda -1}+{\mathit{b}}^{\lambda }\right), \lambda =1, \cdots , \Lambda$$

(45)

Here, the superscript $\lambda$ denotes the layer number of the neural network, ${\mathit{L}}^{\lambda -1}$ the input vector from the $\left(\lambda -1\right)\mathrm{th}$ layer (specially, ${\mathit{L}}^0={\mathit{h}}_t$ and ${\mathit{L}}^{\Lambda }=y_t$), ${\mathit{W}}^{\lambda }$ and ${\mathit{b}}^{\lambda }$ respectively the weight matrices and biases vector of $\lambda \mathrm{th}$ layer. Note that all hidden states of the last GRU are input into the same FNN in turn for training. The total parameters $N^{\mathit{F}}$ in FNN can be calculated as

$$N^{\mathit{F}}=\sum _{\lambda =1}^{\Lambda }{N}^{\lambda }\left(N^{\lambda -1}+1\right)$$

(46)

where $N^{\lambda }$ is the number of neurons of the $\lambda \mathrm{th}$ layer.

3.3. Training Procedure

The loss function, defined as the mean absolute error (MAE), is expressed as

$$\mathcal{l}\left(\mathit{\theta }\right)={\displaystyle \frac{1}{N^D}}\sum _{k=1}^m\left|{\mathit{y}}_k-{\widehat{\mathit{y}}}_k\right|$$

(47)

where ${\mathit{y}}_k$ and ${\widehat{\mathit{y}}}_k$ denote the ground-truth and predicted outputs at the $i\mathrm{th}$ time sequence, $N^D$ is the total number of datapoints in the training dataset, and $\mathit{\theta }$ represents the parameters of the neural networks. These parameters are updated by minimizing the loss function, as follows:

$${\mathit{\theta }}^{\ast }=\underset{\theta }{argmin}\left(\mathcal{l}\left(\mathit{x},\mathit{\theta }\right)\right)$$

(48)

where ${\mathit{\theta }}^{\ast }$ represents the values of $\mathit{\theta }$ that minimize $\mathcal{l}\left(\mathit{x},\mathit{\theta }\right)$. Common optimization methods for minimizing the loss function include gradient decent techniques (such as Adam) and quai-Newton methods (such as L-BFGS). In this work, the Adam optimizer [40] is tentatively employed to update the hyperparameters. The Backpropagation (BP) technique and Backpropagation Through Time (BPTT) technique are used to compute the gradients of the loss function with respect to the parameters of FNN and GRU, respectively. In the following sections, the GRU–FNN model consists of three stacked GRU layers with $100$ units each and a single time-distributed dense layer with $100$ neurons. This configuration balances computational cost and error minimization through trial and error. The leaky rectified linear unit (LeakyReLU) is used as the activation function for the GRU, while the linear function is used for the FNN.

4. Results and Discussions

4.1. Initialization of GRU-FNN Parameters by Training Theoretical Data

The inputs and outputs are unrolled through $100$ time steps, i.e., $\mathit{x}=\left\{{\mathit{x}}_1,{\mathit{x}}_2,···,{\mathit{x}}_m\right\}$ and $\mathit{y}=\left\{{\mathit{y}}_1,{\mathit{y}}_2,···,{\mathit{y}}_m\right\}$ with m = 100. Let $t_k$, ${\lambda }_k$ and $P_k$ denote the loading time, the stretch and the first P-K stress at the $k\mathrm{th}$ time step, respectively. $t_k$ is calculated by the following formula:

$$t_k=\left\{\begin{array}{c}{\displaystyle \frac{{\lambda }_k-1}{\left|\dot{\lambda }\right|}}, \dot{\lambda }>0\\ {\displaystyle \frac{{\lambda }_{ma\mathit{x}}-1}{\left|\dot{\lambda }\right|}}+{\displaystyle \frac{{\lambda }_{ma\mathit{x}}-{\lambda }_k}{\left|\dot{\lambda }\right|}}, \dot{\lambda }<0\end{array}\right.$$

(49)

where ${\lambda }_{ma\mathit{x}}$ is the maximum of $\lambda$ during loading tests. The inputs ${\mathit{x}}_k=\left(t_k,{\lambda }_k\right)$ and the outputs ${\mathit{y}}_k=\left(P_k\right)$ are used for the uniaxial loading-unloading viscoelasticity. Ten samples are generated by capturing 100 points per curve in Figure 2 and then split into training data (80% of total samples) and testing data (20% of total samples). The training data with a batch size of 2 are used during model’s training process, while the testing data are used for validation of the model’s generalization capabilities.

The GRU–FNN model is implemented in the software library Keras 2.10.0 and its information is listed in

Table 1. Information of the GRU-FNN model.

Layer	Output Shape	Parameters
GRU 1	(None, 100,100)	31,200
GRU 2	(None, 100,100)	60,600
GRU 3	(None, 100,100)	60,600
Dense	(None, 100,1)	101

. In the output shape of layers, the first None refers to the batch, the second number the time steps and the third number the units (i.e., the dimension of the hidden state) of GRU or the dimension of the outputs of the dense layer. According to Equation (43), the numbers of parameters of GRU 1 and GRU 2 (GRU 3) layers are calculated respectively as 3 × 100 × (2 + 2 + 100) = 31,200 and 3 × 100 × (2 + 100 + 100) = 60,600. According to Equation (46), the number of parameters of the dense layer is calculated as 1 × (100 + 1) = 101.

Figure 5 illustrates the evolution of loss with respect to epoch during training and testing. It is observed that the loss decreases dramatically at first and then stabilizes at a steady state. Additionally, the testing loss is slightly higher than the training loss. These observations suggest that the model does not suffer from overfitting or underfitting issues.

Figure 6 compares the stress–stretch curves from the training data with the model’s predictions at different stretching rates after 5000 epochs. It was observed that the loss value of the training data stabilizes around 0.50 after approximately 5000 epochs through trial and error. Therefore, we trained the model for 5000 epochs. To evaluate the model’s accuracy under various conditions, the root mean square error (RSME) is employed, defined as $\mathrm{RMSE}={\displaystyle \frac{1}{m}}{{\displaystyle \sum }}_{k=1}^m\sqrt{{\left({\mathit{y}}_k-{\widehat{\mathit{y}}}_k\right)}^2}$, where ${\mathit{y}}_k$ and ${\widehat{\mathit{y}}}_k$ represents the ground truth and the predicted values. The RMSE quantifies the discrepancy between the predicted results and the training or testing data. RMSE values for all cases in this work are listed in Table 2. For the stretching rates of $\left|\dot{\lambda }\right|$ = 0.01/s, 0.02/s, 0.06/s, 0.08/s, 0.12/s, 0.14/s, and 0.16/s, the corresponding RMSE values are 0.76 0.56, 0.53, 0.71, 0.85, 0.83, 0.73, and 0.63, respectively. The predictions from the GRU–FNN model are highly consistent with training data, and the hysteresis curves exhibit typical viscoelasticity.

Figure 7 compares the stress–stretch curves from the testing data with the model’s predictions at two different stretching rates after 5000 epochs. For the stretching rates of $\left|\dot{\lambda }\right|$ = 0.04/s and 0.18/s, the corresponding RMSE values are 3.41 and 1.52, respectively. The testing data, which was not used during the training process, is accurately predicted by the trained model, demonstrating the GRU–FNN model’s strong generalization capabilities in predicting viscoelastic behaviors. The parameter values (weights and biases) here after 5000 epochs are saved for use as initial parameters in training the GRU–FNN model on experimental data.

4.2. Initialization of GRU–FNN Parameters by Training Theoretical Data of Viscoelasticity of VHB4905 with Scarce Data from Experiments

In the thermo-viscoelastic experiments of Liao and Hossain [41], VHB4905 samples with dimensions of 100 m × 10 mm × 0.5 mm are used for cyclic loading–unloading tests to investigate the time- and stretching rate-dependent viscoelasticity of VHB polymers. Three datasets, selected from the experimental stress–stretch curves at 273 K and three different stretching rates, i.e., 0.03/s, 0.05/s, and 0.10/s, are split into two training datasets for learning the viscoelastic behaviors and one testing dataset for checking the GRU–FNN model’s capabilities in predicting viscoelasticity of materials with scarce training data.

First, the three datasets are used with the GRU–FNN model initialized with default parameter values (weights and biases), and the predicted stress is obtained after 5300 epochs. Figure 8 shows the evolution of the loss with respect to epochs during training and testing. It is observed that the loss for testing data are significantly higher than that for training data after several epochs and fluctuates at a high level, indicating model overfitting. Figure 9 compares the stress–stretch curves from the experimental data with the model’s predictions at 273 K after 5300 epochs. For the stretching rates of $\left|\dot{\lambda }\right|$ = 0.03/s, 0.05/s, and 0.10/s, the corresponding RMSE values are 2.84, 3.07, and 62.49, respectively. The training data from cyclic experiments at $\left|\dot{\lambda }\right|$ = 0.03/s and $\left|\dot{\lambda }\right|$ = 0.05/s are predicted accurately. However, the testing data from cyclic experiments at $\left|\dot{\lambda }\right|$ = 0.10/s show substantial divergence with the predicted results, further demonstrating the poor generalization capabilities of the model in predicting viscoelasticity of materials with scarce training data.

To address the overfitting issues associated with scarce data, the parameter values from training the GRU–FNN model with theoretical data are used as the initialization for training on the experimental data. Figure 10 shows the evolution of the loss with respect to epochs during both training and testing. It is evident that the gap between the loss for training and testing data has significantly narrowed compared to Figure 8. Given that the theoretical data were trained for 5000 epochs, the experimental data were trained 300 epochs to maintain consistency in the total training duration (5300 epochs). Figure 11 compares the stress–stretch curves from the experimental data with the model’s predictions at 273 K after 300 epochs. It was also observed that the loss values of training data and testing data stabilize after approximately 300 epochs through trial and error. For the stretching rates of $\left|\dot{\lambda }\right|$ = 0.03/s, 0.05/s, and 0.10/s, the corresponding RMSE values are $1.90$, $2.62$, and $15.70$, respectively. The high degree of consistency between the predicted and experimental stress–stretch curves indicates that our approach has effectively alleviated the overfitting issues.

Finally, we split four datasets, selected from the stress–stretch curves at four different stretching rates, i.e., 0.025/s, 0.05/s, 0.10/s, and 0.20/s in the experiments of VHB4905 samples with dimensions of 130 mm × 10 mm × 0.5 mm [42], into three training datasets and one testing dataset to further assess the GRU–FNN model’s capabilities in predicting the viscoelasticity of materials across different temperatures. The saved parameter values from training with theoretical data are used for initialization. Figure 12 shows the evolution of the loss with respect to epochs during training and testing at 313 K. It can be seen that the gap of the loss during training and testing is smaller than that observed in Figure 9 due to the inclusion of an additional training dataset. Figure 13 compares the stress–stretch curves from the experimental data and the model’s predictions at 313 K after 300 epochs. For the stretching rates of $\left|\dot{\lambda }\right|$ = 0.025/s, 0.05/s, 0.10/s, and 0.20/s at 313 K, the corresponding RMSE values are 0.81, 1.32, 1.72, and 4.55, respectively. Figure 14 shows the evolution of the loss with respect to epochs during training and testing at 333 K. The evolution trend is similar with that in Figure 12. Figure 15 compares the stress–stretch curves from the experimental data and the model’s predictions at 333 K after 300 epochs. At 333 K, for the same stretching rates, the RMSE values are 0.24, 0.31, 0.53, and 4.27, respectively. The model accurately predicts the training data and provides reasonable predictions for the testing data, demonstrating improved generalization capabilities.

Table 2. Statics of RMSE values.

Cases	Stretching Rates (/s)	RMSE Values (kPa)
Figure 6	(0.01,0.02,0.06,0.08,0.10,0.12,0.14,0.16)	(0.76,0.56,0.53,0.71,0.85,0.83,0.73,0.63)
Figure 7	(0.04,0.18)	(3.41,1.52)
Figure 9	(0.03,0.05,0.10)	(2.84,3.07,62.49)
Figure 11	(0.03,0.05,0.10)	(1.90,2.62,15.70)
Figure 13	(0.025,0.05,0.10,0.20)	(0.81,1.32,1.72,4.55)
Figure 15	(0.025,0.05,0.10,0.20)	(0.24,0.31,0.53,4.27)

Table 2. Statics of RMSE values.

Cases	Stretching Rates (/s)	RMSE Values (kPa)
Figure 6	(0.01,0.02,0.06,0.08,0.10,0.12,0.14,0.16)	(0.76,0.56,0.53,0.71,0.85,0.83,0.73,0.63)
Figure 7	(0.04,0.18)	(3.41,1.52)
Figure 9	(0.03,0.05,0.10)	(2.84,3.07,62.49)
Figure 11	(0.03,0.05,0.10)	(1.90,2.62,15.70)
Figure 13	(0.025,0.05,0.10,0.20)	(0.81,1.32,1.72,4.55)
Figure 15	(0.025,0.05,0.10,0.20)	(0.24,0.31,0.53,4.27)

4.3. Analyzing the Sensitivity of the GRU-FNN Model

To test the sensitivity of the GRU–FNN model to input data, a random number is generated from a standard normal distribution (mean = 0, standard deviation = 1). The generated noise value is then multiplied by 0.5% or 1.0% of the original stretch values, resulting in a noise adjustment value that is proportional to the magnitude of the original data. This value is added to the original data to simulate noise in the data. Since the time series can be calculated based on the stretching strain and stretching rate, we only apply noise processing to the stretching strain sequence here. The predicted results of the dataset at 273 K with added noise are shown in Figure 16. It can be observed that the model is resilient to small amounts of noise. However, as the noise increases, the prediction error also increases. When the noise exceeds the error between adjacent points, the prediction curve tends to become unsmooth and may even change direction, leading to a complete failure of the prediction. The same pattern was also observed in the predictions of the other two datasets at 313 K and 333 K, respectively, which are not shown here, but the corresponding RMSE values are listed in

Table 3. Statics of RMSE values.

Cases with Noise	Stretching Rates (/s)	RMSE Values (kPa)
273 K_0.5%	(0.03,0.05,0.10)	(1.93,2.17,8.28)
273 K_1.0%	(0.03,0.05,0.10)	(2.16,3.80,20.30)
313 K_0.5%	(0.025,0.05,0.10,0.20)	(0.45,0.59,0.84,4.55)
313 K_1.0%	(0.025,0.05,0.10,0.20)	(1.84,1.23,0.83,4.73)
333 K_0.5%	(0.025,0.05,0.10,0.20)	(1.37,0.84,0.66,5.27)
333 K_1.0%	(0.025,0.05,0.10,0.20)	(0.50,0.54,0.59,5.16)

.

5. Conclusions

In this paper, a physics-guided GRU–FNN model for predicting the viscoelasticity of solids at large deformation is proposed and validated by the comparison of the stress–stretch curves from the model prediction and the uniaxial experiments of VHB4905. The major novelty of the present work lies in the following aspects. First, the time- and loading path- dependent stress can be predicted by utilizing both time and stretch/strain sequences as inputs of the proposed GRU–FNN model. Second, the saved values of the parameters after training theoretical data are employed as the parameter initialization of GRU–FNN model for learning the viscoelastic behavior of solids with scarce experimental data, which significantly alleviates the overfitting issues that arise in the case of scarce data. The results show that the GRU–FNN model can reduce the RMSE by several times, significantly improving its generalization ability in scarce data scenarios. The physics-guided initialization of parameters in this work provides a new idea on how to incorporate physical knowledge into a machine learning model when the governing equations describing physical phenomenon in one material cannot be specified. The model can be applied to predict multiaxial stress–strain relations of viscoelastic solids in a similar manner when the experimental data from multiaxial tests are available, which will be our future study subject. This study facilitates the rapid prediction of materials’ nonlinear viscoelastic properties and lays a theoretical foundation for future real-time monitoring and performance prediction of polymer-based devices such as electronic artificial muscles and dielectric elastomer actuators.