## 1 Introduction  
With the advancement of materials science, numerous polymer materials have emerged in engineering [1; 2], including rubber, plastics, nylon, organic glass, textile fibers, and other industrial materials. Additionally, geological materials such as asphalt, concrete, soil, rocks, and biological materials like bones and muscles fall into this category. These materials exhibit the deformation characteristics of elastic solids and the time-dependent flow characteristics of viscous fluids during operation and use, collectively referred to as viscoelastic materials [3]. The mechanical behavior of such materials under loading depends on the magnitude of the load and varies with time. Typically, they can be regarded as linear viscoelastic materials under small strain levels. When conducting finite element analysis on linear viscoelastic materials, there are two distinct implementation approaches based on the chosen form of constitutive relationship. One approach involves employing a differential form of constitutive relationship [4]. Thisapproach represents the viscoelastic mechanical behavior through various combinations of force-bearing elements, such as springs and dashpots. Different combinations yield different constitutive models [5], such as the generalized Maxwell model [6], the generalized Kelvin model [7; 8], and so forth. This constitutive relationship is expressed in the following form of a differential expression:  
\[P[\sigma(t)]=Q[\varepsilon(t)] \tag{1}\]  
where \(P\) and \(Q\) represent linear differential operators with respect to time \(t\), and \(P=P(D)=\sum\limits_{r=0}^{n}p_{r}\frac{\dot{\sigma}}{d\dot{r}^{\prime}}\), \(Q=Q(D)=\sum\limits_{r=0}^{n}q_{r}\frac{\dot{\sigma}}{d\dot{r}^{\prime}}\). The expanded form of the differential constitutive relationship is denoted as \(\sigma+p_{1}\dot{\sigma}+p_{2}\ddot{\sigma}+\cdots=q_{0}\varepsilon+q_{1} \dot{\varepsilon}+q_{2}\ddot{\varepsilon}+\cdots\), where \(p_{r}\) and \(q_{r}\) are material constants obtained through experiments. In different constitutive models, specific coefficients in the equation may be zero. The existing literature often employs a combination of a differential constitutive relationship and the initial strain method in developing viscoelastic finite element analysis programs, typically using the generalized Kelvin model. This model consists of a spring element and a series of Kelvin bodies (each comprising a spring and a dashpot connected in parallel). This form readily separates the total strain at any given moment into elastic strain (provided by the spring element) and viscous strain (the sum of strains from each Kelvin body). Consequently, the viscous strain can be regarded as the initial strain, enabling the use of the initial strain method employed in elastic finite element analysis for the solution. Specifically, time is discretized into a series of intervals, assuming that stress remains constant within each interval. Based on the characteristics of the generalized Kelvin model, the stress on individual spring elements and each Kelvin body connected equals the known stress within that interval. Using [PERSON]'s creep compliance [9], the viscous strain of each Kelvin body under this stress is calculated and summed to obtain the total viscous strain, and this total viscous strain is treated as the initial strain. Additional forces are determined, and displacement increments caused by viscous strain during that interval are calculated. The overall analysis results are calculated step-by-step. This method maintains the stiffness matrix unchanged and equal to the initial stiffness matrix when solving the linear equation set at each interval, thus exhibiting higher computational efficiency. However, using the initial strain method mentioned above relies on two underlying conditions: firstly, the ability to express the total strain at any given moment as the sum of elastic strain and viscous strain; secondly, the assumption that the material's Poisson's ratio remains constant.