Finite Element Analysis for Linear Viscoelastic Materials Considering Time-Dependent Poisson’s Ratio: Variable Stiﬀness Method

: For linear viscoelastic materials, this paper proposes a ﬁnite 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 stiﬀness matrix is established, thereby forming a convolutional ﬁnite element equation. Furthermore, the recursive calculation of the convolutional integral is derived, and the calculation steps of the ﬁnite 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 veriﬁed through two numerical example s with linear viscoelastic material. The results indicate that the proposed variable stiﬀness method for the ﬁ nite 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 stiﬀness method proposed in this paper eﬀectively avoids the assumption of constant stress during the micro time interval, thus signiﬁcantly enhancing the ﬁnite 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.


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 Citation: Wang, X.; Gao, J.; Wang, Y.; Bai, J.; Zhao, Z.; Luo, C. Finite Element Analysis for Linear Viscoelastic Materials Considering Time-Dependent Poisson's Ratio: relationship.One approach involves employing a differential form of constitutive relationship [4].This approach 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: where P and Q represent linear differential operators with respect to time t , and , where r p and r q 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 Kelvin'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.The first condition is only satisfied when employing structures similar to the generalized Kelvin model.This implies the necessity of knowing the creep compliance beforehand, which requires providing corresponding data through creep tests.The second condition typically applies only to some incompressible materials.However, meeting both conditions in practical engineering can be challenging due to many problems.For example, in the viscoelastic analysis of solid rocket motor propellants [10,11], which are polymer materials, relaxation tests rather than creep tests are usually conducted.These tests provide data on the relaxation modulus varying with time instead of the creep compliance [12].Additionally, Poisson's ratio is not constant.In such cases, the direct application of the initial strain method is not feasible.
Another approach is to employ an integral constitutive relationship [13].This approach considers the stress-time history sum of a series of rectangular pulses over small time intervals.Utilizing the Boltzmann superposition principle, the following form of the hereditary integral is derived [14]: where ( ) σ t and ( ) ε τ represent stress and strain, respectively, both of which are func- tions of time; ( ) E t represents the relaxation modulus of the viscoelastic material; and ε 0 represents the initial strain.Current commercial software packages such as ABAQUS 2022 and ANSYS 2023 utilize this constitutive relationship [15,16].However, these commercial programs require the input of shear relaxation modulus ( ) G t and bulk relaxation modulus

( )
K t , which decompose deformation at any given time into volumetric and shear components.Typically, these two relaxation moduli are assumed to follow a Prony series, with material parameters fitted from experimental data.Once these relaxation modulus expressions are determined, the remaining process is similar to the elastic finite element method.Time is discretized into a series of intervals, and numerical integration methods are used to obtain displacement and strain for each time interval.The advantage of this method lies in its generality.Regardless of the specific constitutive model used, as long as the Prony series of the relaxation modulus is obtained, the analysis can be completed following steps similar to those of traditional finite element methods.However, the drawback is apparent: forming the stiffness matrix must be repeated for each time step, resulting in high computational complexity and lower efficiency.Another drawback is the requirement for shear and bulk relaxation tests to obtain corresponding material parameters for these commercial software packages.Since these tests are operationally complex, they are rarely conducted in practical engineering.Most relaxation tests involve tension or compression [17], providing the relaxation modulus ( ) E t over time [18,19].Then, assuming a constant Poisson's ratio, the bulk and shear moduli are derived from these parameters and input into the commercial software.However, this conversion process inevitably introduces additional errors.Therefore, while this method uses shear relaxation and bulk relaxation moduli, it relies solely on the elastic relaxation modulus and a constant Poisson's ratio [20][21][22].
In summary, it is clear that whether developing custom finite element programs or using commercial software for analysis, it is commonly assumed that Poisson's ratio remains constant throughout the entire analysis process.However, this assumption is inaccurate for many practical viscoelastic problems.For example, experimental results on Poisson's ratio of viscoelastic materials like polymers show that it is not constant throughout the testing.It tends to increase towards 1/2 with rising temperature and prolonged loading time, while it decreases towards 1/3 with decreasing temperature and shorter loading time.Even slight variations in Poisson's ratio can significantly impact the calculation results of mechanical properties [23,24].Therefore, this paper aims to develop a viscoelastic finite element algorithm that can simultaneously account for the relaxation behavior of the elastic modulus and the creep Poisson's ratio.Compared with the initial strain method, the calculation accuracy is improved by considering the time variation in Poisson's ratio.Compared with commonly used commercial software, the required constitutive parameters are easier to obtain, thus providing a foundational algorithm for addressing such practical engineering problems.

The Generalized Maxwell Model for Viscoelastic Materials
Without loss of generality, the generalized Maxwell model is used to represent the relaxation properties of the elastic modulus in viscoelastic materials.Defining the elastic relaxation modulus [25][26][27]: The creep properties of viscoelastic materials are defined using the generalized Kelvin model, and the creep Poisson's ratio is established [28,29]: where E ∞ represents the elastic coefficient of the spring element in the generalized Max-

Integral Constitutive Relationship of Viscoelastic Materials
Based on the elastic relaxation modulus and creep Poisson's ratio previously mentioned, the elastic relaxation matrix ( )  can be defined as follows: With the elastic relaxation matrix above, the integral constitutive relationship of viscoelastic materials [30] can be expressed as follows: where D t represents the elastic relaxation matrix at time t , { } ε 0 represents the in- itial strain vector, and ( ) { } ε τ represents the strain vector at time τ .After obtaining the constitutive relationship of viscoelastic materials, similar to the finite element derivation process, the relationship between the nodal forces and nodal displacements of the element can be derived using the principle of virtual displacement: ( ) ( ) the relaxation stiffness matrix of the element at time t , where     B represents the strain matrix.
Integrating these individual matrices yields the global relaxation stiffness matrix after obtaining the relaxation stiffness matrix for any element at a specific time.Simultaneously, by assembling the elemental load vectors into the global load vector and introducing the boundary conditions, the global equation at that time can be derived as follows: where ( )

{ }
P t represents the global load vector at time t .

The Foundation of Convolutional Integral for Viscoelastic Materials
Numerical implementation is required to solve Equation ( 8) involving convolutional integrals.Here, the recursive format for convolutional integration is given.Suppose that we have the following form of a convolutional integral equation, where 1 t represents the initial moment of integration and t represents the end moment of integration.In practical scenarios, If the integration time t t , then we have j small time intervals, thus According to the mean value theorem for integrals, the summation can be simplified as follows [31]: Further simplification yields the following: Substituting the above equation into the original convolutional integral Equation (9) yields the following: After rearranging the above equation, the following recursive equation for solving the convolutional integral equation is derived: . Subsequently, starting from the third time point ( ≥ 2 j ), +1 j y 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 ① Discretize the time into an array as follows: where     t j represents the j -th element in the array, with ② Determine material parameters: use experimental data and the least squares to fit the material parameters in Equations ( 3) and ( 4): ③ Determine the corresponding load vector for each time point: ( ② Use the relaxation matrix obtained above and derive the initial stiffness matrix for each element at the initial time: ③ Assemble the global stiffness matrix . It is important to note that the stiffness matrix at the initial time must be saved for subsequent calculations.
④ Establish the global load vector at the initial time ⑥ Solve the following equations to obtain the displacement { } (3) Calculation at the second time point, i.e., =   =   1 2 t t t and ( ) and further construct the relaxation matrix for this time point: ③ Use the relaxation matrix obtained above to formulate the stiffness matrix for each element: ④ Assemble the global stiffness matrix ( ) ⑥ Establish the global load vector for this time point  (4) From the third time point onwards, repeat the following steps: ④ Assemble the global stiffness matrices ( ) ⑥ Calculate the additional force using the following steps: Initialize the additional force vector { } Formulate the relaxation matrices Calculate the stiffness matrices for each element.
Formulate the global stiffness matrices end ⑦ Establish the load vector ( ) { } P t j     at the nodes for this time point.
at this time point, i.e., =     t t j , and save this displacement.

Analytical Solution
The analytical solution can be obtained by solving the general expression of the differential constitutive equation of linear viscoelastic materials as follows: ( ) ( ) where ( ) Node 2 where 1 E and 1 η are the elastic and viscous coefficients of the Maxwell elements, re- spectively; E ∞ is the elastic coefficient of the spring, representing the elastic modulus of the material as time approaches infinity.
Referring to the general expression (16), for the H|M model, we have the following: ( ) ( ) Rewriting Equation ( 16) yields the following: ( ) ( ) It can be seen that ( ) ( ) According to the principle of elasticity-viscoelasticity correspondence, the solutions to linear viscoelastic problems are exactly the same as those of elastic problems under identical boundary conditions in Laplace space.Simply substitute E in the elastic solution with the Laplace space expression ( ) ( ) Here,

( )
Q s is the Laplace transform of ( ) Q D , and ( ) P s is the Laplace transform of ( ) P D .For the H|M model used in the given example ( ) ( ) Therefore, ( ) ( ) For the axial tension rod, its elastic displacement solution is . According to the principle of elasticity-viscoelasticity correspondence, the linear viscoelastic displacement solution retains the same form in Laplace space.To obtain the viscoelastic displacement solution in Laplace space, simply substitute E with ( ) ( ) in Equation ( 21), thus resulting in the expression shown below: ( ) Organize Equation ( 22) to obtain a form that is easy to perform inverse Laplace transform with as follows: After this arrangement, each term in Equation ( 23) can be identified using the Laplace transform table to find its corresponding inverse transformation form.The analytical solution of the viscoelastic problem in the time domain is then derived by applying the inverse Laplace transform to Equation ( 23), as shown below: It can be organized as follows: ( 1 , In the given example, 1 10.8 MPa.s η = ; in the case of 1 3.6 s τ = , we have the following: (

Comparison and Discussion of Simulation Results
Following the calculation steps previously outlined, considering the case of only one element, the viscoelastic displacements for different material parameters can be obtained.As depicted in Figure 2, the graph also includes the analytical solution obtained using Equation (25).As presented in Figure 3, the relative error of two methods under different material parameter values can be calculated using the following formula: where numeric u is the axial tension displacement obtained by the numerical algorithm proposed in this study, analytical u is the analytical solution to the displacement obtained according to Equation (25), and error E is the relative error of the results obtained by the two methods.From Figure 2, it is evident that in the case of using only one element, the results yielded by the algorithm proposed in this paper nearly perfectly align with the analytical solution.Under constant loading, the viscoelastic displacement of the axial rod progressively increases with time, indicating the creep properties inherent to viscoelastic materials.Analyzing the error variation depicted in Figure 3, it becomes apparent that initially, the error is minimal, followed by an increase over time, reaching a peak around the relaxation time.Subsequently, the error gradually diminishes.Throughout the process, the error consistently remains negative, indicating a slight underestimation in the calculated results compared to the analytical solution.The initial ascent and subsequent descent of the error can be attributed to the rapid change in the slope of the curve near the relaxation time, as evident from the equation for the relaxation modulus of the generalized Maxwell model.This suggests a steep slope of the curve in proximity to the relaxation time, succeeded by a more gradual change beyond this range, signifying a relatively stable relaxation modulus.In this example, for a relaxation time of 3.6 s, the stress decreases by 63.2% from its initial value ( 1 e of the initial stress), indicative of rapid attenuation.Conse- quently, the material parameters undergo rapid changes near the relaxation time (or within a time range of the same order of magnitude as the relaxation time).Due to the uniform time intervals used in the calculation, the error increases within this interval.However, beyond this range, as the elastic modulus changes more slowly, the error decreases.To enhance the calculation accuracy, reducing the time intervals near the relaxation time and enlarging them away from the relaxation time can be contemplated.

Information Detail
As depicted in Figure 4, we consider a straight rod with an equal cross-section characterized by a unit weight.(i.e., the corners of the bottom free end) is calculated in Figure 4.

Analytical Solution of Viscoelastic Displacement
According to elastic mechanics, the analytical solution for the vertical displacement of elastic materials is as follows: ( ) where E and µ are the elastic modulus and Poisson's ratio of the material, respectively.Using the same approach as in Section 4.1 and based on the principle of elasticityviscoelasticity correspondence, the solution of Equation (28) in the phase space after Laplace transform can be obtained as follows: and ( ) ( ) Based on this framework, we derive the following: Substituting the above equation into Equation (30) and introducing 0 1 yields the following: ( ) ( ) Applying the inverse Laplace transform to Equation (35) yields the following: Similarly, applying the inverse Laplace transform to Equation (31) yields the following: ( ) ( ) From this, the inverse Laplace transform of Equation ( 29) can be derived as follows: ( ) where ( ) ( ) In this example, the analytical for viscoelastic displacement, accounting for both Poisson's ratio and elastic modulus variations over time, can be obtained by substituting the respective parameters as follows:

= − = − ×
Based on the algorithm proposed in this paper, the viscoelastic displacement of the sectioned straight bar under self-weight is calculated.The time interval is set to 1 s, and the finite element mesh is partitioned, as shown in Figure 5.The calculation results of the viscoelastic displacement at the free end corner within 50 s are shown in Figure 6, where the analytical solution obtained by Equation ( 41) is also provided.From the calculation results, it can be observed that for the sectioned straight bar under self-weight, at the position of the free end corner, the instantaneous elastic displacement due to loading is approximately 0.5 mm.Subsequently, due to the viscous flow of the material, the viscoelastic displacement gradually increases with time, reaching 2.84 mm at 50 s, which is approximately five times the initial displacement.A comparison with the analytical solution shows that the maximum error of the algorithm proposed in this paper is about 0.3%, indicating good computational accuracy.From the above two examples, it can be observed that when deriving analytical solutions for viscoelastic problems while leveraging the principle of elasticity-viscoelasticity correspondence, it is necessary to first obtain the analytical solution for the corresponding elastic problem.Subsequently, by applying Laplace transform, the analytical solution for the viscoelastic problem in the Laplace domain can be derived.Following this, by performing an inverse transformation, the analytical formula for the viscoelastic problem can be obtained.In many practical scenarios, however, obtaining analytical solutions for corresponding elastic problems is challenging due to the complexity of boundary conditions, geometric shapes, and loading conditions involved.Consequently, the applicability of analytical solution methods based on Laplace transform is severely limited.In contrast, the algorithm proposed in this study falls under the finite element method, which is not constrained by the conditions mentioned above.Therefore, it proves to be more suitable for addressing complex engineering problems encountered in practice.

Conclusions
This paper introduces a viscoelastic finite element algorithm that simultaneously considers the relaxation behavior of elastic modulus and the creep Poisson's ratio in viscoelastic materials.The relaxation behavior of the elastic modulus is modeled using the generalized Maxwell model, while the creep Poisson's ratio is represented using the generalized Kelvin model.Starting from the integral constitutive relationship and incorporating the principle of virtual displacement, a convolution-based variable stiffness method for solving equations is proposed.The paper also presents recursive calculations for convolution-based numerical integration.Based on this, the calculation steps of the proposed algorithm are outlined.The main conclusions of this paper can be summarized as follows: (1) In contrast to the initial strain method derived from the differential constitutive relationship, this algorithm can simultaneously account for the variation in the elastic modulus and Poisson's ratio over time, thus establishing a more precise simulation pathway and a novel perspective for modeling the mechanical properties of viscoelastic materials.(2) Compared to the initial strain method for linear viscoelastic materials, the proposed variable stiffness method considering the time-dependent Poisson's ratio effectively avoids the assumption of constant stress over small time intervals, thereby significantly enhancing the accuracy of the finite element calculations for linear viscoelastic materials.(3) Compared to the finite element method commonly used in commercial software, the proposed variable stiffness method considering the time-dependent Poisson's ratio, with the elastic modulus and Poisson's ratio as material parameters, establishes a simulation algorithm that aligns with experiments on viscoelastic materials.

)
Elastic calculation at the initial time ( = 0 t) relaxation matrix at the initial time according to the following equation:

⑤
Introduce any necessary boundary conditions.
point, i.e., the second time point, and save this displacement value.

Discussion of the Example 4 . 1 .
Viscoelastic Displacement of the Rod Subject to Uniaxial Tension 4.1.1.Parameter Specification A vertical axial rod is fixed at the bottom and subjected to tensile force P , with ma- terial behavior described by the generalized Maxwell viscoelastic model, comprising a spring element and a Maxwell body in parallel.It is assumed that only one-dimensional elements are used.The load 2 .8M Pa s ⋅ , and 21.6 M Pa s ⋅ .The corresponding relaxation times are 1 1.8 τ = , 3.6 and 7.2 s, as shown in Figure 1.The visco-elastic displacement of the rod subject to the uniaxial tension at different relaxation times is solved in this section, where the numerical results are computed using our own programmed algorithms.

Figure 1 .
Figure 1.Axial rod of the linear viscoelastic material.
generalized Maxwell model, composed of a spring element and a Maxwell element in parallel, degenerates into the traditional Poynting-Thomson model, commonly known as the H|M model (H represents Hookean spring, M represents a Maxwell element, and H|M represents a spring and a Maxwell element in parallel).Its differential constitutive equation is as follows:

Figure 2 .
Figure 2. Comparison of different relaxation times.

Figure 3 .
Figure 3. Relative error at different relaxation times.

.
The rod is fixed at the upper end and free at the lower end, with a length of L = 2 m (equivalent to 2000 mm) and a cross-section of 400 mm 400 mm × .The origin is denoted as O, with the z-axis extending positively upward and the y-axis oriented horizontally and positively to the right.Considering the material as linear viscoelastic and using the generalized Maxwell model to reflect the relaxation characteristics of the material's elastic modulus as elastic modulus is 0 3.65MPa E =.Using the generalized Kelvin model to reflect the creep Poisson's ratio of the material, where 0 s ratio at the final stable state is 0.49 µ ∞ =.The relaxation time is 3.6 s, i.e., = = .The description of the material pa- rameters is explained in the previous example.

Figure 4 .
Figure 4. Slender members with equal cross-section and viscoelastic material.The vertical displacement under the self-weight effect at points located at 0 z = and 200 mm x y = = the Laplace transform of Poisson's ratio elastic modulus in phase space, respectively.For this example, within the generalized Maxwell model describing the relaxation characteristics of elastic modulus, m = 1, and in the generalized Kelvin model used to describe the creep Poisson's ratio, n = 1.Consequently, we have the following:

Figure 5 .
Figure 5. Infinite element mesh for the cross-section of straight bars.

Figure 6 .
Figure 6.Displacement of a straight bar with equal cross-section under self-weight effect.
Kelvin body in the generalized Kelvin model, and τ i represents the relaxation time of the i -th Kelvin body in the generalized Kelvin model.m and n represent the number of Maxwell bodies in the generalized Maxwell model and the number of Kelvin bodies in the generalized Kelvin model, respectively.The parameters above are fitted based on data from the relaxation test for the elastic modulus and the creep test for Poisson's ratio.
represents the global relaxation stiffness matrix at time t ; { } δ t represents the global displacement vector at time t ; and j , and ( )