For linear viscoelastic materials, this paper proposes a finite element analysis method based on an integral constitutive relationship that can simultaneously consider the relaxation behavior of the elastic modulus and the creep Poisson’s ratio. Firstly, the generalized Maxwell model is employed to depict the relaxation characteristics of the elastic modulus, while the generalized Kelvin model is used to represent the creep Poisson’s ratio. Subsequently, the element relaxation stiffness matrix is established, thereby forming a convolutional finite element equation. Furthermore, the recursive calculation of the convolutional integral is derived, and the calculation steps of the finite element for viscoelasticity considering the time-dependent nature of both the elastic modulus and Poisson’s ratio are established. Finally, the accuracy of the proposed algorithm is verified through two numerical examples with linear viscoelastic material. The results indicate that the proposed variable stiffness method for the finite element analysis of linear viscoelastic materials can simultaneously consider the changes in the elastic modulus and Poisson’s ratio over time, thereby establishing a new path and idea for the more accurate simulation of viscoelastic materials’ mechanical properties. Compared with the initial strain method for linear viscoelastic materials, the variable stiffness method proposed in this paper effectively avoids the assumption of constant stress during the micro time interval, thus significantly enhancing the finite element calculation accuracy of linear viscoelastic materials. The proposed method establishes a simulation algorithm that matches existing commercial software with viscoelastic material experiments by considering the elastic modulus and Poisson’s ratio as material parameters.