## 1 Introduction
### Integral Constitutive Relationship of Viscoelastic Materials  
Based on the elastic relaxation modulus and creep Poisson's ratio previously mentioned, the elastic relaxation matrix \([D(t)]\) can be defined as follows:  
\[[D(t)]=\frac{E(t)}{(1+\mu(t))(1-2\mu(t))}\begin{bmatrix}1-\mu(t)&\mu(t)&\mu(t)& 0&0&0\\ \mu(t)&1-\mu(t)&\mu(t)&0&0&0\\ \mu(t  
equation: \(y_{2}=y_{1}+\frac{2\mu(t_{2})}{A(0)+A(t_{2}-t_{1})}\). Subsequently, starting from the third time point (\(j\geq 2\)), \(y_{j+1}\) can be solved. Calculation Steps of Finite Element Analysis for Viscoelastic Materials Considering Time-Dependent Poisson's Ratio  
Based on the recursive calculation provided by Equation (15) to obtain the convolutional integral, the viscoelastic finite element Equation (8) given in Section 2 can be solved. After rearrangement, the calculation steps for the variable stiffness method, which considers both the relaxation behavior of the elastic modulus and the creep Poisson's ratio, are as follows:  
1. Data Initialization 2. Discretize the time into an array as follows: \(t= \begin{bmatrix}0,&t_{1}&t_{2}&t_{3}&\cdots&t_{N-1}\end{bmatrix} _{1\times N^{ \prime}}\) where \(t[j]\) represents the \(j\)-th element in the array, with \(t[1]=0\), \(t[2]=t_{1}\), \(t[3]=t_{2}\), \(t[j+1]=t_{j}\); additionally, \(\Delta t_{j}=t[j+1]-t[j]=t_{j}-t_{j-1}\) is introduced. 3. Determine material parameters: use experimental data and the least squares to fit the material parameters in Equations (3) and (4): \(m\), \(E_{\infty}\), \(E_{r}\), \(\tau_{r}\) (\(r=1,2,\cdots,m\)), \(n\), \(\mu_{0}\), \(\mu_{i}\), \(\tau_{i}\) (\(i=1,2,\cdots,n\)). 4. Determine the corresponding load vector for each time point: \(P=\begin{\{}P_{0}&P_{1}&\cdots&P_{N-1}\end{\}}\). 2. Elastic calculation at the initial time (\(t=0\)) 3. Calculate \(E(0)=E_{0}=E_{\infty}+\sum\limits_{r=1}^{m}E_{r}\) and \(\mu(0)=\mu_{0}\) and derive the relaxation matrix at the initial time according to the following equation: \[[D(0)]=\frac{E_{0}}{(1+\mu_{0})(1-2\mu_{0})} \begin{bmatrix}1-\mu_{0}&\mu_{0}& \mu_{0}&0&0&0\\ \mu_{0}&1-\mu_{0}&\mu_{0}&0&0&0\\ \mu_{0}&\mu_{0}&1-\mu_{0}&0&0&0\\ 0&0&0&\frac{1-2\mu_{0}}{2}&0&0\\ 0&0&0&0&\frac{1-2\mu_{0}}{2}\end{bmatrix} \]  
Use the relaxation matrix obtained above and derive the initial stiffness matrix for each element at the initial time: \([k(0)]^{e}=\int\limits_{v}{[B]^{T}[D(0)][B]dv}=\int\limits_{-1}^{1}\int \limits_{-1}^{1}\int\limits_{-1}^{1}{[B]^{T}[D(0)][B]\cdot|J|d\delta d\eta d\zeta}\). 4. Assemble the global stiffness matrix \([K(0)]\). It is important to note that the stiffness matrix at the initial time must be saved for subsequent calculations. 5. Establish the global load vector at the initial time \(\{P(0)\}\), denoted as \(P_{0}\). 6. Introduce any necessary boundary conditions. 7. Solve the following equations to obtain the displacement \(\{\delta\}_{1}\) at the initial time \(t[1]\). 8. \(\left[K(0)\right]\cdot\left\{\delta\right\}_{1}=\{P(0)\}\) 9. Calculation at the second time point, i.e., \(t=t[2]=t_{1}\) 10. Calculate \(\Delta t_{1}=t[2]-t[1]=t_{1}-0=t_{1}=t[2]\). 10. Calculate \(E(\Delta t_{1})=E_{\infty}+\sum\limits_{r=1}^{m}E_{r}e^{-\frac{\Delta t_{1}} {\tau_{r}}}\) and \(\mu(\Delta t_{1})=\mu_{0}+\sum\limits_{i=1}^{n}\mu_{i}\bigg{(}1-e^{-\frac{ \Delta t_{1}}{\tau_{i}}}\bigg{)}\) and further construct the relaxation matrix for this time point: \[\begin{array}{llll}[D(\Delta t_{1})]&=\frac{E(\Delta t_{1})}{(1+\mu( \Delta t_{1}))(1-2\mu(\Delta t_{1}))}\\ \times& \begin{bmatrix}1-\mu(\Delta t_{1})&\mu(\Delta t_{1})&\mu(\Delta t_{1} )&0&0&0\\ \mu(\Delta t_{1})&1-\mu(\Delta t_{1})&\mu(\Delta t_{1})&0&0&0\\ \mu(\Delta t_{1})&\mu(\Delta t_{1})&1-\mu(\Delta t_{1})&0&0&0\\ 0&0&0&\frac{1-2\mu(\Delta t_{1})}{2}&0&0\\ 0&0&0&0&\frac{1-2\mu(\Delta t_{1})}{2}\end{bmatrix}