On relevance of time-dependent Poisson’s ratio for determination of relaxation function parameters

The current study concerns the determination of material constants of a three-dimensional linear viscoelastic model. It is assumed that the constitutive equation utilizes a Prony series as a memory function. A method for the evaluation of relaxation function parameters is presented which can be used for arbitrary loading histories. The proposed methodology is applied to the identification of the viscoelastic constants of acrylonitrile butadiene styrene (ABS). For that purpose, a number of rheological tests in tension have been performed on ABS standard dogbone specimens. The significance of the time-dependent Poisson’s ratio for the determination of material parameters is investigated. It is found that taking into account the measurements of specimen’s lateral contraction over time has a particularly strong influence on the identified values of parameters responsible for the bulk behavior. Several boundary value problems have been analyzed in order to assess the influence of the material parameter values on the obtained solutions. It is demonstrated that some oversimplifications assumed during the determination of viscoelastic constants can lead to a loss of precision or even wrong results.


Introduction
The Prony series is probably the most commonly used generic function of viscoelasticity. From the mechanistic point of view, this function yields from the generalized Maxwell rheological model [10,13,29,37]. The constitutive equations of linear viscoelasticity utilizing the Prony series are implemented in numerous finite element analysis programs, such as Abaqus [17], ADINA [1] or MSC Marc [25]. This is the reason why the evaluation of the viscoelastic constants of this particular model becomes a compelling issue.
In practice, it is common to examine the properties of a viscoelastic solid using the creep tests, e.g., [24]. This approach is justified by the fact that standard creep experiments require relatively simple apparatus compared to other testing techniques. Another reason is that a uniaxial stress state during a creep test results in a simplification of the one-dimensional process equations which yield from the compliance form of tensor constitutive equation, e.g., [37]. However, the widely used engineering software is mainly based on the displacement formulation of the finite element method (FEM). This fact requires the compliance constitutive equation of viscoelasticity to be inverted in order to obtain the parameters of the relaxation function which could be further utilized in the stiffness constitutive equations used by the finite element software.
In the case of material memory functions using more than one exponential term, the inversion of the constitutive equation utilizing an analytical approach is usually troublesome, if possible. Thus, numerous numerical algorithms have been proposed to enable interconverting generic functions which employ multiple exponential terms. Such algorithms are usually based on the relation between creep and relaxation functions which can be expressed in the form of a convolution integral, e.g., [3,12]. A different approach is to utilize Technical Editor: João Marciano Laredo dos Reis. the iterative adjusting of the measured and theoretical complex moduli [2,21].
Through the years, numerous algorithms for the determination of the Prony series material parameters were proposed. Due to the presence of the exponential terms, the generic functions of viscoelasticity are highly sensitive to the parameter values. Garbarski [14,15] utilized a systematic search algorithm to determine the constants of double exponential function from multistep rheological test data for various thermoplastic polymers. Chambers [6] and Chen [7] independently employed pattern search algorithms to evaluate the constants of multiple exponential functions from relaxation data including ramp loading. The more rapid, yet less stable gradient methods are occasionally used for the identification of viscoelastic constants, e.g., [3].
Some researchers postulate that better convergence and stability are achieved if the material characteristic times are excluded from the classical optimization. This can be achieved by setting the number of relaxation times and their values a priori, e.g., [6]. On the other hand, Ciambella et al. [9] introduced a method where the relaxation times were estimated using a quasi-genetic algorithm, whereas the viscoelastic coefficients were iteratively optimized using both gradient and pattern search methods. A modification of the aforementioned algorithm was proposed by Suchocki [32]. It appears that currently, due to the increased computational power of modern computers, the relaxation times in both linear and nonlinear viscoelasticity do not require any special treatment and can be determined by means of the least squares method the same way other material constants are evaluated. This fact allows for a reduction in number of the utilized exponential terms.
The material parameter evaluation algorithms discussed in the literature concern in major part one-dimensional rheological models; thus, many problems occurring during the identification of three-dimensional, tensor constitutive models are skipped. It should be emphasized that in the case of relaxation data, the material constants determined for a one-dimensional prototype of a constitutive equation are usually in no way usable for the tensor equations utilized by FEM. This fact is associated with the existence of three-dimensional strain state and a time-dependent Poisson's ratio which have to be taken into account during the derivation of the process equations. The aforementioned problem was discussed by Qvale and Ravi-Chandar [30] who proposed approximating the measured lateral strain data with a polynomial function of time in order to enable evaluating the hereditary integral and obtaining a closed-form process equation. This relationship describing the stress relaxation could be further used for the purpose of approximating the relaxation data.
An analytical integration of the hereditary integrals usually results in lengthy and complicated process equations which have to be further handled during the determination of material parameters. What is more, the closed-form result of evaluating a hereditary integral can be obtained only for relatively simple deformation histories. Goh et al. [16] demonstrated for a nonlinear viscoelastic model that exact values of the material constants can be identified when a numerical integration of the constitutive equation is used instead of a closed-form process equation. For that purpose, the numerical integration algorithm developed by Taylor et al. [33] was utilized. The same procedure of material parameter evaluation was later successfully applied by Sorvari and Malinen [31] in the case of a one-dimensional linear viscoelastic model.
In this work, a newly developed method for determining the material constants of a three-dimensional linear viscoelastic model is introduced. The material memory function is assumed in the form of a Prony series. The presented identification algorithm utilizes the Taylor's numerical integration method to calculate the theoretical material response which is further used during the least squares optimization procedure. The algorithm can be used for arbitrary strain histories. The nonconstant Poisson's ratio which generates the effect of lateral strain evolution during a uniaxial relaxation test can be taken into account easier than it was proposed before [30,37]. Thus, the parameters of the shear and bulk relaxation functions can be evaluated directly from the data collected in a standard uniaxial stress relaxation experiment. There is no need for performing any other more elaborate rheological test, subsequent identification of the material parameters and further conversion of the determined function of viscoelasticity in order to obtain the parameters of the shear and bulk relaxation functions. The proposed material parameter identification algorithm is applied to evaluate the viscoelastic constants of acrylonitrile butadiene styrene (ABS) on which the uniaxial stress relaxation tests in tension have been performed. The collected measurements include the specimen's lateral contraction changing over time which results in a timevariable strain deviator and axiator. The significance of the lateral strain's measurement for the determination of viscoelastic constants is investigated. A number of boundary value problems are analyzed in order to assess the influence of the identified material parameter values on the obtained solutions.

Basic notions
The stiffness form of the constitutive equation of isotropic linear viscoelasticity takes the form [10,13,37]: with tr (•) being the trace operator. The shear and bulk relaxation functions are assumed in the form of Prony series, i.e.
where G ∞ is the long-term shear modulus, G j and j ( j = 1, 2, … , N ) are the coefficients and relaxation times, respectively, accounting for the nonequilibrium shear response, whereas B ∞ is the long-term bulk modulus, B j and ̃j ( j = 1, 2, … , M ) are the coefficients and relaxation times, respectively, accounting for the nonequilibrium bulk response.
Sometimes it is convenient to rewrite Eq. (2) in an equivalent differential form by utilizing the so-called viscoelastic shear and bulk nonequilibrium stresses. i.e., j ( j = 1, 2, … , N ) and h j ( j = 1, 2, … , M ), respectively. Thus: Equations (5) 2 and (5) 4 can be derived based on the analysis of mechanistic Maxwell model, e.g., [10,13,37], whereas the additivity of nonequilibrium stresses yields from the generalized Maxwell model (Fig. 1) which may be also called the H-nM 1 model. By taking advantage of the Laplace transform, i.e. (3) Eq. (2) can be written in a form which resembles the stress-strain relationships of the linear elasticity: The fictitious "elastic constants" Ḡ � (s) and B � (s) can be introduced in the form of the following functions of the transformation parameter s: which are called the shear and bulk relaxances, respectively [35], and are related to the so-called stretch relaxance Ē � (s) and Poisson's retardance ̄�(s) by the following formulas: with where ̄a(s) and ̄a(s) are the Laplace transforms of the axial stress and axial strain, respectively, whereas ̄l(s) is the transform of the lateral strain in an infinitesimal cube subjected to uniaxial tension. The time-dependent Young's modulus and Poisson's ratio are given by the following equations: which yields In the case of a strain excitation defined by the Heaviside's step function H(t), i.e., a (t) = 0 H(t) , the time-dependent elastic constants are given as [18,35]: The relationships given above facilitate the method of solving the boundary value problems called the correspondence principle, e.g., [10,13,37]. This method assumes that when a solution to a problem of the elasticity theory is known, a solution to a corresponding problem of the viscoelasticity theory can be found by replacing the elastic constants by the proper relaxances and retardances and subsequent inverse Laplace transform. The applicability of this method is limited to the cases when the boundary conditions have the form of a multiplication of independent time and coordinate functions.

Model computation algorithm
The numerical algorithm which is utilized for the time integration of Eq. (2) is a modification of the numerical scheme introduced by Taylor et al. [33]. The calculations performed during a single time step are listed in the table below. An assumed strain history is treated as an input for the algorithm. The shear and bulk relaxation functions are taken in the form of Prony series as given in Eq. (4). The presented algorithm is utilized to compute the theoretical deviatoric and volumetric stresses which are used for the calculation of the total square error during the determination of material parameter values [16,32].

Numerical integration of linear viscoelastic model
Calculate strain deviator and dilatation from time step n + 1 3. Calculate total Cauchy stress in time step n + 1 4. Store stresses and strains:

Discretized scalar equations for uniaxial stress relaxation
In order to determine the scalar process equations, it is assumed that during the stress relaxation an infinitesimal material cube is subjected to a uniaxial stress state (i.e., 11 = 11 (t) , 22 = 33 = 0 ) and a three-dimensional strain state. The components of the strain tensor in a Cartesian coordinate system Ox 1 x 2 x 3 can be gathered in a square matrix: which have to be solved in order to determine the theoretical stress components 11 (t) and s 11 (t) along with the volumetric stress p(t). It should be emphasized that while 0 is an imposed kinematic excitation, the lateral strain l (t) is an unknown time function which has to determined experimentally. Qvale and Ravi-Chandar [30] used a polynomial function to approximate the lateral strain measurements. By this way, an analytical function l (t) was obtained which, after inserting into Eqs. (17) 1 and (17) 3 , allowed to solve the hereditary integrals.
In common practice, a simplification based on the assumption of a constant Poisson's ratio = const is utilized, i.e. Substitution of Eq. (18) into Eqs. (15) and (16)  It follows from Eqs. (17) 2 , (17) 4 and (20), that in the case of constant Poisson's ratio, the shear and bulk relaxation functions differ by a constant factor only, i.e., However, postulating = const a priori without any experimental justification can lead to a loss of precision or even serious errors, as it will be demonstrated in the sections to follow. The problems of solving the hereditary integrals (17) 1  and (17) 3 can be skipped when the constitutive equation of linear viscoelasticity is integrated numerically by utilizing the algorithm presented in the previous paragraph. For the input time histories of axial and lateral strains, the axial component of deviatoric strain and the dilatation at the time step t n+1 are given as: while the deviatoric viscoelastic overstress component is defined by the following recurrence-update formula: whereas the volumetric overstress The total axial component of the stress deviator at time step t n+1 is given as: while the volumetric stress Equations (21)(22)(23)(24)(25) are utilized within the material parameter determination algorithm to compute the theoretical stress response for the experimentally measured histories of e 11 (t) and (t) . The experimental strain histories have been obtained by performing the uniaxial stress relaxation tests on acrylonitrile butadiene styrene (ABS).

Experimental setup
For all the conducted mechanical tests of ABS, the material testing machine Zwick/Roell Z005 was used with the maximum load capacity of 5 kN. The samples were gripped into the pressure jaws with constant pressure of about 4 MPa (Fig. 2). For the purpose of assessing the real elongation of polyamide specimens, the long-range extensometer with the total measurement range of about 1.2 m was employed.
In order to simultaneously measure both the specimen's elongation and its lateral contraction, an optical method based on the digital images correlation (DIC) was utilized. The method was invented in the 1980s [4,8,28] and has a wide range of applications in the mechanical testing of materials [5,22,23]. The DIC method allows to obtain the maps of displacements and macroscopic deformations measured on a surface of a tested specimen. A DIC stand consists of a monochromatic digital camera with a suitable lens and a PC with software (Fig. 3a). The images of the specimen's surface were collected by the DIC system at the frequency of 1 Hz. A commercial program Vic2D [36] was later used for the postprocessing of images. Before the test, random spots with high contrast must be applied to the surface of the specimen (Fig. 3b, c). After selecting the field for analysis, it is divided into smaller subareas characterized by a unique distribution of shades of gray. The computer algorithm compares digital photographs of a nondeformed sample with the photographs taken after the deformation and finds the new positions of analyzed subareas by the method of best fit. As a result, a map of displacements or macroscopic deformations can be obtained for the entire analyzed area. The accuracy of displacement measurements depends on the camera's resolution and reaches 0.01 pixels.
The dogbone specimens with a rectangular cross section ( 4 × 10 mm) and the reduced section length of 60 mm were used for the experiments. The samples were designed according to the guidelines of PN-EN ISO 527-2:2012 The initial phase of the research is comprised of several introductory relaxation tests which were conducted in order to check the performance of the DIC method and verify that the results obtained for the same rheological experiments are repeatable. Additionally, uniaxial ramp tension experiment was performed for the purpose of measuring the basic mechanical characteristics of the examined material. A yield stress y = 40.36 MPa was registered for the axial strain a = 2.51% . The specimen's failure was observed for the stress f = 45.82 MPa and the axial strain max = 43.74% . These measurements are in a good agreement with the specification provided by the producer [34]. The strain field at the gauge length of the specimen remained homogenous until failure and no necking were observed.
The data collected during the uniaxial ramp test were utilized to specify the parameters of the main relaxation experiment. For the purpose of mimicking the case of a Heaviside step excitation, a preloading was performed at the strain rate set to ̇= 0.01 s −1 with the target value of the axial strain set to 2% . After the maximum axial strain of 2.56% was reached, the testing machine's automatic control system corrected the strain value. Subsequently, the axial strain of 2% was held constant for 3600 s allowing the axial stress to relax. The straining period lasted 4.5 s which is short compared to the total stress relaxation time. This fact justifies the assumption that the measured material response can be viewed as a response to a Heaviside step loading. In order to eliminate any initial transient effects which might have been introduced during preloading, the data from the initial 8 s of the experiment were discarded from the analysis.

Determination of material parameters
The experimental measurements collected during the stress relaxation test are presented in Fig. 4, i.e., the axial strain 11 (t) = a (t) , the lateral strain 22 (t) = l (t) and the axial stress 11 (t) . Equations (21) were utilized to calculate the time history of the axial strain deviator component e 11 (t) and the dilatation (t) which are shown in Fig. 5. The collected data were used for the determination of viscoelastic constants of ABS.
The material parameters were evaluated by utilizing the least squares method. For that purpose, the "fminsearch" function available in Scilab software was utilized. The following total square error function was defined for the minimization: is the column matrix of the searched parameters and s 11 ( ) and s 11 are the theoretical and experimental axial deviatoric stress, respectively, whereas p( ) and p are the theoretical and experimental volumetric stress, respectively, at the time instants k ( k = 1, 2, … , K ). The collocation points were chosen to be distributed uniformly in the logarithmic time. The long-term elastic moduli G ∞ and B ∞ were calculated from the following relations: with the approximations 11 (∞) ≈ 11 (t max ) and Two different approaches were employed for the identification of viscoelastic constants. In the first approach, the experimental time histories of the axial strain deviator component and the dilatation, measured during the relaxation test (Fig. 5), were used as an input for Eqs. (22)(23)(24)(25) to compute the theoretical values of stresses in Eq. (26). In this case, the time-dependent Poisson's ratio was taken into account. Two versions of the viscoelastic model were considered, i.e., the simplified version ( N = 1 and M = 0 in Eq. (4) which leads to the standard solid model in shear) and the extended version ( N = 2 and M = 2).
The comparison of the experimental results and the simplified model predictions ( N = 1 , M = 0 ) can be seen in Fig. 6 while the parameter values are gathered in Table 1 (model 1).   In the second approach, it was assumed that the material parameters of viscoelasticity would be determined based on the measurements of the axial strain and the axial stress, solely (simulated case when the lateral strain data are unavailable because a proper measurement was not performed). Thus, the specimen's lateral contraction is taken to be a linear function of the axial strain, i.e., l (t) = − a (t) , where a (t) = 0 H(t) while is a constant value. This assumption corresponds to the case of a constant Poisson's ratio. Based on the information which can be found in the available technical documentation, it is assumed that = 0.35 [11,34] The comparison of the experimental results and the theoretical stresses according to the simplified model ( N = 1 , M = 0 ) is shown in Fig. 8 while the parameter values are listed in Table 1 (model 2). The curve fitting achieved using the extended model ( N = 2 , M = 2 ) can be viewed in Fig. 9 with the parameter values gathered in Tables 2 and 3 (model 2).
Since the loading is defined by the Heaviside step function, the time-dependent Young's modulus and Poisson's ratio can be found using the so-called quasi-elastic approximation [29], i.e.

Equations (28) and (4) were utilized to calculate E(t) and
(t) for the material parameter values gathered in Tables 2  and 3. The constant sets labeled as "Model 1" were used for that purpose. The experimentally measured time-dependent Young's modulus and Poisson's ratio were computed using Eq. (13). It can be seen in Fig. 10 that a very good agreement was found between the experimental data and the functions calculated from Eqs. (28) and (4). Both the experimentally measured and the simulated Poisson's ratios  are decreasing with time. The Poisson's ratio of polymers is generally reported as an increasing function of time, e.g., [27]. The observed decrease is possibly caused by the secondary relaxations which are associated with the motions of side groups or short segments of the main polymer chains.

Analysis of selected boundary value problems
In order to investigate the influence of the material parameter identification method on the obtained results, several examples are analyzed including two variations of the Lamé problem. Below both the analytical and the approximate solutions are discussed.

Thick-walled tube
As an example of a two-dimensional problem, let us consider a thick-walled tube with the external radius of a = 100 mm and the internal radius of b = 50 mm (Fig. 11). The tube is loaded with the internal pressure p = 23 MPa; thus, the boundary conditions in the polar coordinates are [13]: where rr is the radial stress component. The elastic solution of the considered problem is given by the following relationships: and with tt being the tangential stress component while u r and u t are the radial and the tangential displacement components, respectively. It should be noted that according to Eq. (30), the stresses are independent of the assumed material properties.
In the case of simplified viscoelastic model which assumes a single Prony term for simulating the shear response and an elastic volumetric behavior ( N = 1 , M = 0 ), an analytical viscoelastic solution can be found by utilizing the correspondence theorem. For that purpose, the stress-strain relations given by Eq. (2) have to be written in the form with p ′ 1 , q ′ 0 , q ′ 1 and q ′′ 0 being the functions of material parameters. After applying the Laplace transform according to Eq. (6), it is found that: . (39) The change of the radial displacement over time as predicted by the analytical solution given by Eqs. (41) and (42) is plotted in Fig. 12 for the two different material parameter sets (model 1 and model 2) which are gathered in Table 1.
Since the pressure loading is defined by the Heaviside step function, an approximate solution to the analyzed problem can be found. The so-called quasi-elastic approximation [29] is found by replacing elastic constants in Eq. (31) 1 by the proper time functions. After substituting Eq. (28) into Eq. (31) 1 the following approximate solution is found: The radial displacement's time evolution according to the quasi-elastic approximation is plotted in Fig. 12. It can be seen that for both material parameter sets gathered in Table 1, the analytical solution and the quasi-elastic approximation are in a very good agreement.
In the case of extended viscoelastic model for ABS ( N = 2 and M = 2 in Eq. (4)), obtaining an analytical solution encounters a problem since the inverse Laplace transform of Eq. (40) 1 could not be found analytically. The quasielastic solutions for both material parameter sets collected in Tables 2 and 3 are plotted in Fig. 13.

Pin press fitted into infinite plate
A variant of the Lamé boundary value problem was considered where a rigid pin is press fitted into a hole in infinite linear viscoelastic plate (Fig. 14). The radius of the hole is r 0 = 50 mm while the radius of the pin is greater by the value of u 0 = 1.5 mm. The problem can be solved using polar coordinates. In the case of quasi-elastic approximation, the radial and tangential stresses are given by the following equations [26]: whereas the displacement field is given as . and is independent of the material properties of the plate. The radial and tangent stresses given by Eq. (44) are plotted in Fig. 15. The shear relaxation function utilizing two Prony terms was assumed along with the material parameter values gathered in Table 2.

Cyclic volumetric deformation of cube
In the third of the analyzed problems, a linear viscoelastic cube was considered which is subjected to a cyclic dilatation. The strain measures are given as follows: The dilatation (t) is defined by a trapezoidal function which is shown in Fig. 16a. Since the components of the strain deviator are constant and equal to zero, the material response  is purely volumetric. The theoretical predictions of the cyclic volumetric stress p(t) for the two different material parameter sets (Table 3) are shown in Fig. 16b.

Conclusions
In this work, a method for the determination of material parameters of the three-dimensional linear viscoelasticity has been presented. The method allows to evaluate the viscoelastic constants for an arbitrary deformation history. In the case of a uniaxial stress relaxation test, the proposed methodology enables taking into account the true stress deviator and dilatation histories as depicted in Fig. 5. This in turn results in a more precise values of the identified material parameters than in the case of making a common assumption about the constant Poisson's ratio.
What is more, it has been observed that a simplification based on assuming the relaxation times a priori, e.g., [6,7] or applying any special systematic search techniques for the purpose of determining them [9,32] is unnecessary at least for the case of linear viscoelasticity. The computational power of the modern computers allows to identify the values of the relaxation times along with the other constitutive constants by taking advantage of the available optimization packages such as Scilab or MATLAB, for instance.
The discussed material parameter identification algorithm was applied to determine the viscoelastic constants of ABS thermoplastic polymer. Two different linear viscoelastic models were analyzed. The simplified model assumes a linear elastic bulk behavior and a single Prony term for describing the shear response. The extended model for ABS utilizes two Prony terms for capturing both the shear and the bulk behavior (independently). The simplified version of the constitutive model allows for an acceptable description of the material's shear response (Figs. 6a and 8a). Some analytical solutions of the boundary value problems can be obtained for this model. The extended constitutive model is able to simulate the material response in both shear and bulk deformations with a very good precision (Figs. 7 and 9). In the case of ABS polymer, adding additional Prony terms to the extended version of the model proved pointless since no improvement in the curve fitting was achieved this way.
Determining the material parameters of linear viscoelasticity with the assumption of a constant Poisson's ratio does not strongly affect the identified values of the relaxation times, cf Tables 1, 2 and 3. However, the errors in the evaluated values of the bulk relaxation function coefficients are substantial, cf Table 3. This fact significantly affects the solutions of problems involving volumetric deformations (Fig. 16b).
The obtained results lead to a conclusion that utilizing material parameters which were determined assuming a constant Poisson's ratio does not necessarily have to generate completely incorrect solutions to the boundary value problems, cf Figs. 12, 13 and 15. Nevertheless, this simplification will always cause a certain loss in precision. Generally, the uniaxial stress relaxation data used for the determination of the viscoelastic constants should include the measurements of the specimen's lateral contraction, if possible. It is particularly important for the correctness of solutions to boundary value problems where a significant volume change takes place.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.