Numerical and Experimental Investigations of Polymer Viscoelastic Materials Obtained by 3D Printing

The aim of this research is to determine the relaxation and creep modulus of 3D printed materials, and the numerical research is based on the finite volume method. The basic material for determining these characteristics is ABS (acrylonitrile butadiene styrene) plastic as one of the most widely used polymeric materials in 3D printing. The experimental method for determining the relaxation functions involved the use of a creep test, in which a constant increase of the stress of the material was performed over time to a certain predetermined value. In addition to this test, DMA (dynamic mechanical analysis) analysis was used. Determination of unknown parameters of relaxation functions in analytical form was performed on the basis of the expression for the storage modulus in the frequency domain. The influence of temperature on the values of the relaxation modulus is considered through the determination of the shift factor. Shift factor is determined on the basis of a series of tests of the relaxation function at different constant temperatures. The shift factor is presented in the form of the WLF (Williams-Landel-Ferry) equation. After obtaining such experimentally determined viscoelastic characteristics with analytical expressions for relaxation modulus and shift factors, numerical analysis can be performed. For this numerical analysis, a mathematical model with an incremental approach was used, as developed in earlier works although with a certain modification. In the experimental analysis, the analytical expression for relaxation modulus in the form of the Prony series is used, and since it is the sum of exponential functions, this enables the derivation of a recursive algorithm for stress calculation. Numerical analysis was performed on several test cases and the results were compared with the results of the experiment and available analytical solutions. A good agreement was obtained between the results of the numerical simulation and the results of the experiment and analytical solutions.


Introduction
Most polymeric materials are used today for their good mechanical properties in relation to manufacturing costs. For this reason, the mechanical properties of these materials are of interest because of their wide application. It is important to know the dependence of deformation and stress when analyzing the stability of elements or products exposed to a load, and it is known that the characteristics of polymeric materials change over time. A complete description of the mechanical properties of polymeric materials usually requires knowledge of the behavior of the material over 10 to 20 years [1]. This means that in order to determine the reliability and durability of a certain polymer structure, it is necessary to know the mechanical characteristics during the entire time of use.
In the theory of linear viscoelasticity, such characteristics are most often the relaxation modulus E(t) and the creep function J(t), which define the response of the material to the applied load (i.e., the Heaviside function). Moreover, due to the short duration of the experiment, it is often necessary to replace the step load function with other non-standard modules from the frequency domain to time domain was proposed by Schapery [3], while Emri and Tschoegl [24] developed a recursive algorithm that generates a line spectrum from experimental data. They applied the window method to ensure positive values of all elements in the Prony series. Ramkumar used the regularization method with quadratic programming to minimize oscillations in the experimental data. Later, Schapery and Park developed a new method of interconversion based on a numerical and approximateanalytical approach. To solve the problem of determining the coefficients in the Prony series, many tools are available in the literature, and some of these have already been implemented in commercial software packages such as MATLAB, Mathcad, Maple, Origin, etc. [1][2][3]. In this paper, the MATLAB software package with built-in subroutines for optimization was used.
One of the most important methods for determining the long-term behavior of polymers is based on the principle of temperature-time superposition (TTSP), which uses information from short-term tests at different temperatures. The method uses the similarity of the variations of the relaxation/creep modulus over time at different temperatures. This approach allows the creation of a master curve that moves based on the value of the shift factor. Later, the experimental time frame was extended to successfully show the relationship between time, creep/relaxation function values, and glass region temperature. One-dimensional tensile experiments were performed using a Meissner-type rheometer to measure ABS properties and a material model was proposed that takes into account the effects of strain rate, solidification strain, and temperature sensitivity. An innovative thermoforming process has recently been developed to form composite layers using pressure while melting the polymer during the injection phase. This process is an integration of thermoforming and injection processes, and can result in a significant reduction in the tools, machines, and operating costs needed. A very significant study of the mechanical properties of 3D printed polymeric materials was conducted by a group of researchers [25], whose subject was initially wood as a natural polymer, to extend the research to other synthetic materials (PLA, ABS, PE) [26][27][28] as well as their combinations.
Based on the results of previous research, it is possible to introduce an analytical expression in the form of a Prony series to express the values of the relaxation function during an arbitrary period of time. The unknown parameters in the Prony series can be determined by regression analysis, and the influence of temperature can be represented using the temperature-time principle of superposition. Analytical expressions derived in this way can be used in numerical analysis, and the obtained results are confirmed on the basis of conducted experiments or available analytical solution. The model developed by Demirdzic, Dzaferovic and Ivankovic [13] for the analysis of thermoviscoelastic deformations in the finite volume method is the basis for calculations, while a recursive algorithm obtained on the basis of using Prony's series is an upgrade that enables more efficient and faster calculation of parameters.

Experimental Analysis
This chapter describes the experimental analysis of viscoelastic materials obtained by 3D printing. The material characteristics and basic tests for determining relaxation functions are described, as well as the tests used in this paper.
Tests with a constant increase in load over time can be considered a limit case when the time of the ramp function coincides with the total load time. Here, a constant rate of change of load over time is used instead of a constant load. Tests with a constant increase in load are typical for testing soft materials, such as gels or liquids. The most interesting and practical uses are optimization methods where relaxation functions are obtained from the process of approximating the stress or deformation of the material. Relaxation functions are represented as a sum of exponential functions, i.e., a Prony series [24]. The appropriate mathematical expression for the Prony series is: where E ∞ , E k , τ k are the "rubbery" relaxation modulus, coefficients in the Prony series, and relaxation times, respectively. While the appropriate formulation for the creep modulus is: where J o , J k , and τ k are the coefficients in the Prony series and relaxation times in case of creep. On the other hand, the constitutive relation for a linear viscoelastic material is given in the following form: Based on Equation (1), the following expressions can be obtained for the relaxation modulus: or the creep compliance σ o are determined from the approximation or nonlinear regression procedure described in the next section.
In cases when we have the experimentally determined value of the relaxation module in tabular and/or graphical form, the task is to determine the analytical form of the relaxation module E rel (t). The most efficient and common method for approximating this type of result is the least squares method or its variation. This problem boils down to determining the parameters E k and τ k in Equations (4) and (5). If the difference of the squares between the measured values and the calculated values is expressed in the following form: where L is the number of measurements. Then, from the condition that the stated sum is minimum, we have a system of nonlinear equations: From the last expression (7) arises a system of nonlinear equations with as many unknowns as we have parameters in the Prony series (E k , i τ k ) whose solution gives an analytical expression for the relaxation module that can be further used in numerical calculations. Numerous nonlinear optimization routines are already available in the MATLAB Optimization Toolbox. In this paper, the "lsqcurvefit" MATLAB functions were used to solve the problem of nonlinear approximation (fitting), which in turn use the least squares method or its variation with an interface designed specifically for data approximation problems.
Relaxation and creep tests are suitable for the long-term study of material responses (minutes or days), but less accurate in cases of shorter times (seconds or less). Dynamic tests are used when the stress (or strain) resulting from sinusoidal stress is often good for filling a "short" range of polymer reactions. When the viscoelastic material is subjected to a sinusoidally variable stress, the response is also sinusoidal and has the same angular frequency, but is late in the phase with an angle δ. The deformation lags behind the stress at the phase angle δ, and this also applies if the deformation is a controlled variable (input load) instead of the stress. Using Fourier transforms, the real part of the relaxation module can be expressed as follows [1,2]: Then, the imaginary part is: where is k is the number of members in the Prony series. The relaxation modulus of polymeric materials usually has three regions, which are determined by time and temperature.
The first area is called glassy, where the properties are constant and the material is hard and brittle. The test material then enters a transition area where the material properties change drastically. By further increasing the temperature we reach the rubber region, where the material behaves like rubber. The relaxation modulus curve usually extends over several decades, which is physically impossible to measure. For this reason, data from experiments are collected within a certain time interval, which is usually called an experimental window ( Figure 1). frequency, but is late in the phase with an angle δ. The deformation lags behind the stress at the phase angle δ, and this also applies if the deformation is a controlled variable (input load) instead of the stress. Using Fourier transforms, the real part of the relaxation module can be expressed as follows [1,2]: Then, the imaginary part is: where is k is the number of members in the Prony series. The relaxation modulus of polymeric materials usually has three regions, which are determined by time and temperature.
The first area is called glassy, where the properties are constant and the material is hard and brittle. The test material then enters a transition area where the material properties change drastically. By further increasing the temperature we reach the rubber region, where the material behaves like rubber. The relaxation modulus curve usually extends over several decades, which is physically impossible to measure. For this reason, data from experiments are collected within a certain time interval, which is usually called an experimental window ( Figure 1). In this paper, a simple CFS methodology of displacement in closed form is used, which completely eliminates the previously mentioned problems related to the displacement of experimental data in the process of construction of the main curve at the selected reference temperature [3]. Suppose that the experiments were performed at M different temperatures, or {Tj; j = 1, 2, …, M}. Furthermore, let each of the M measured segments be given with Nj discrete data, then it could be written in the following relation: Usually, the experimental data of the relaxation modulus and the creep function are presented as discrete data pairs in a double logarithmic scale. Taking this into account, then the whole range of experimental data can be expressed as: log E , = log E T , log t , = log t T ; j = 1, 2 , 3, … , M; i = 1, 2, 3, … . , N . Figure 1. Schematic representation of a relaxation module for a polymer [24].
In this paper, a simple CFS methodology of displacement in closed form is used, which completely eliminates the previously mentioned problems related to the displacement of experimental data in the process of construction of the main curve at the selected reference temperature [3]. Suppose that the experiments were performed at M different temperatures, or {T j ; j = 1, 2, . . . , M}. Furthermore, let each of the M measured segments be given with N j discrete data, then it could be written in the following relation: Usually, the experimental data of the relaxation modulus and the creep function are presented as discrete data pairs in a double logarithmic scale. Taking this into account, then the whole range of experimental data can be expressed as: E(T, t) = log E j,i = log E i T j , log t j,i = log t i T j ; j = 1, 2 , 3, . . . , M; i = 1, 2, 3, . . . ., N j . If we select the temperature T κ as the reference temperature, then for the corresponding data measured at this temperature as the reference segment of the main curve, the following holds: To construct a master curve at the selected temperature Tκ, it is necessary to move all segments measured at T j > T κ to the right, and the remaining segments, measured at T j < T κ to the left along the logarithmic time axis, so as to obtain a "smooth" main curve, as shown in the schematic in Figure 2. The individual horizontal displacements, log a (T j , j = κ) obtained in this process, represent a set of discrete values describing the so-called time-temperature superposition, which is usually modeled with the WLF expression.
To construct a master curve at the selected temperature Tκ, it is necessary to move all segments measured at Tj > Tκ to the right, and the remaining segments, measured at Tj < Tκ to the left along the logarithmic time axis, so as to obtain a "smooth" main curve, as shown in the schematic in Figure 2. The individual horizontal displacements, log a (Tj, j ≠ κ) obtained in this process, represent a set of discrete values describing the so-called timetemperature superposition, which is usually modeled with the WLF expression.
C1, C2-WLF constants and are determined experimentally, in which To is the reference temperature for which the master curve is constructed, and Ti is the current temperature.
The research samples were printed from Z-ABS ABS filament (Zortrax, Poland) with a Zortrax M200 3D printer. A 0.4 mm nozzle and a 1.75 mm diameter thread were used to print 50 × 12 × 4 mm samples for rheometer tests and 80 × 10 × 4 mm for bending tests. The print nozzle temperature was 275 °C and the printer bay temperature was 80 °C. Samples with different layer thicknesses of 0.09 mm, 0.19 mm, and 0.39 mm were printed. The samples had the same structure-the three outer layers were filled with material, and the inner structure was a mesh structure with a square size of 1.25 mm and a line thickness of 0.4 mm. The rheometer ( Figure 3) "ARES-G2 TA instruments" was used in the research for dynamic mechanical analysis of samples. The specimens ( Figure 4) were attached to the torsional rectangular fixture and the torsional sinusoidal load. The test procedure was controlled by computer Trios software. The amplitude test with oscillating amplitude was performed at 1 Hz with a stress increase of 1.0 × 10−6 to the end of the linear viscoelastic region of the material, which was determined by a sudden change (more than 5%) of the storage modulus. The tests were performed at a constant temperature of 30 °C. The storage modulus, loss modulus, and tan δ are calculated from the material response measurements. C 1 , C 2 -WLF constants and are determined experimentally, in which T o is the reference temperature for which the master curve is constructed, and T i is the current temperature.
The research samples were printed from Z-ABS ABS filament (Zortrax, Poland) with a Zortrax M200 3D printer. A 0.4 mm nozzle and a 1.75 mm diameter thread were used to print 50 × 12 × 4 mm samples for rheometer tests and 80 × 10 × 4 mm for bending tests. The print nozzle temperature was 275 • C and the printer bay temperature was 80 • C. Samples with different layer thicknesses of 0.09 mm, 0.19 mm, and 0.39 mm were printed. The samples had the same structure-the three outer layers were filled with material, and the inner structure was a mesh structure with a square size of 1.25 mm and a line thickness of 0.4 mm.
The rheometer ( Figure 3) "ARES-G2 TA instruments" was used in the research for dynamic mechanical analysis of samples. The specimens ( Figure 4) were attached to the torsional rectangular fixture and the torsional sinusoidal load. The test procedure was controlled by computer Trios software. The amplitude test with oscillating amplitude was performed at 1 Hz with a stress increase of 1.0 × 10 −6 to the end of the linear viscoelastic region of the material, which was determined by a sudden change (more than 5%) of the storage modulus. The tests were performed at a constant temperature of 30 • C. The storage modulus, loss modulus, and tan δ are calculated from the material response measurements.   Oscillatory rheometry is one of the DMA-like techniques used to investigate the mechanical properties of the materials. The material is oscillating, and the response is measured by different methods. The response can be elastic (represented by elastic or storage modulus (G rolled cous (represented by modulus of loss (G″). For real materials the response is a combination of elastic and viscous response. The value tan (δ) is the ratio between modulus of loss and storage module and is sometimes called the dumping factor. Tensile tests of the samples were performed on a universal testing machine Zwick-Roell Z005 with a constant load speed of 1 mm/min.

Numerical Method
This chapter presents a mathematical model to describe the deformations of a viscoelastic solid, including the equation of momentum, the equation of thermal energy, the constitutive relations as well as the initial and boundary conditions for the solution domain [13].

Governing Equations
The law of conservation of momentum and the law of conservation of energy of a viscoelastic solid, of volume V located inside a closed surface S, can be expressed mathematically in the following form [13]:  Oscillatory rheometry is one of the DMA-like techniques used to investigate the mechanical properties of the materials. The material is oscillating, and the response is measured by different methods. The response can be elastic (represented by elastic or storage modulus (G rolled cous (represented by modulus of loss (G″). For real materials the response is a combination of elastic and viscous response. The value tan (δ) is the ratio between modulus of loss and storage module and is sometimes called the dumping factor. Tensile tests of the samples were performed on a universal testing machine Zwick-Roell Z005 with a constant load speed of 1 mm/min.

Numerical Method
This chapter presents a mathematical model to describe the deformations of a viscoelastic solid, including the equation of momentum, the equation of thermal energy, the constitutive relations as well as the initial and boundary conditions for the solution domain [13].

Governing Equations
The law of conservation of momentum and the law of conservation of energy of a viscoelastic solid, of volume V located inside a closed surface S, can be expressed mathematically in the following form [13]: Oscillatory rheometry is one of the DMA-like techniques used to investigate the mechanical properties of the materials. The material is oscillating, and the response is measured by different methods. The response can be elastic (represented by elastic or storage modulus (G rolled cous (represented by modulus of loss (G"). For real materials the response is a combination of elastic and viscous response. The value tan (δ) is the ratio between modulus of loss and storage module and is sometimes called the dumping factor. Tensile tests of the samples were performed on a universal testing machine Zwick-Roell Z005 with a constant load speed of 1 mm/min.

Numerical Method
This chapter presents a mathematical model to describe the deformations of a viscoelastic solid, including the equation of momentum, the equation of thermal energy, the constitutive relations as well as the initial and boundary conditions for the solution domain [13].

Governing Equations
The law of conservation of momentum and the law of conservation of energy of a viscoelastic solid, of volume V located inside a closed surface S, can be expressed mathematically in the following form [13]: where t is the time, ρ is the density, u is the displacement, δσ is the stress increment, n is the outer normal of surface S, δf b is the increment of force, e is the internal energy, q is the heat flux, and s e the source term.

Constitutive Relations
For an isotropic viscoelastic material, the stress history including thermal effects is described by a hereditary integral [13]: where is the volume relaxation modulus, α(T) is the thermal expansion, I is the unit tensor, and T is the body temperature. It is assumed that the relations of the Lame relaxation functions and parameters for elastic materials also apply to viscoelastic materials, and that they are defined as follows: where E(t) is the Young's relaxation modulus, and ν(t) is the Poisson number.
In the case of non-isothermal changes in the viscoelastic material, then times t and τ are replaced by the reduced times using the following next equation: The time temperature shift factor a T (T) is the material characteristic that gives the relationship between the properties of the material at temperature T and the reference temperature T 0 . The relation between the heat flux and the temperature gradient is given by Fourier's law: where k is the thermal conductivity. The relationship between thermal energy and temperature is given by the following expression: where c is the specific heat.

Initial and Boundary Conditions
The equation for a solid is of the hyperbolic type, and at time t = t o the displacement increment and the time gradient of the displacement increment at all points of the solution domain are given as follows [13]: The boundary conditions must be specified on all boundary areas of the resolution domain. They can be of the Dirichlet and/or Von-Neumann type depending on whether the dependent variables or gradient values of the dependent variables are given at the domain boundary.
Dirichlet boundary conditions: where ϕ is T, du. Neumann boundary conditions: In many practical cases, the problem can be solved only in one part of the solution domain, because there is a symmetry of the solution in the remaining area on the axis of symmetry.

Correction of Boundary Conditions
For the case when the increment of force is given as a boundary condition at the boundary of the solution domain, due to the relaxation of the total stress, which follows from the form of the constitutive relation in which relaxation functions are used which are monotonically decreasing, it is necessary to correct the boundary condition as follows [13]: Figure 5 shows the stress relaxation in the case of uniaxial load, when a constant stress σ o is given at the limit. If f 2 (t m ) is the total force at time tm, and df 2 (t m ) is the increment of the force in the m−th time and/or incremental step ( Figure 5), then, using the expression for total stress [15], it holds that: where: the residual increment of force, s B is the boundary surface vector at the domain boundary. Stress σ(t m ) is total stress at time t m−1 which is relaxed at time t m .
Polymers 2021, 13, x FOR PEER REVIEW 9 of 24 Dirichlet boundary conditions: where φ is T, du. Neumann boundary conditions: In many practical cases, the problem can be solved only in one part of the solution domain, because there is a symmetry of the solution in the remaining area on the axis of symmetry.

Correction of Boundary Conditions
For the case when the increment of force is given as a boundary condition at the boundary of the solution domain, due to the relaxation of the total stress, which follows from the form of the constitutive relation in which relaxation functions are used which are monotonically decreasing, it is necessary to correct the boundary condition as follows [13]: Figure 5 shows the stress relaxation in the case of uniaxial load, when a constant stress σo is given at the limit. If f2(tm) is the total force at time tm, and df2(tm) is the increment of the force in the m−th time and/or incremental step ( Figure 5), then, using the expression for total stress [15], it holds that: where: the residual increment of force, s is the boundary surface vector at the domain boundary. Stress σ(tm) is total stress at time tm−1 which is relaxed at time tm.

Calculation of Total Stress
The values of stress tensor σ(tm) are determined in the centers of the control volumes adjacent to the boundary surfaces of the solution domain. The calculation of the total stress can be done using the constitutive relation (16). By introducing the relaxation module from expression (1) into the constitutive relation, we have the following equation:

Calculation of Total Stress
The values of stress tensor σ(t m ) are determined in the centers of the control volumes adjacent to the boundary surfaces of the solution domain. The calculation of the total stress can be done using the constitutive relation (16). By introducing the relaxation module from expression (1) into the constitutive relation, we have the following equation: The first term on the right side of the equation can be considered elastic while the second term is viscoelastic and can be further separated: Individual terms in Equation (27) can be written as: ∂τ dτ is the viscoelastic part of the total stress from the previous time step t m−1 . The remaining members are calculated according to the following relation: The formula for calculating the total stress can be written as follows: It can be seen from relation (29) that to calculate the value of the total stress at the current moment t m it is sufficient to know the values of the viscoelastic part of the stress from the previous time step (σ m−1 visco ) as well as the deformation values from these two time steps.

Experimental Results
This chapter presents the results of processing experiments in which a number of tests were applied, such as the test with constant increase in load (stress/strain), creep test, dynamic frequency test, dynamic amplitude test, and test for determining the temperature shift factor. All these tests were performed on the previously described experimental line and devices that give sufficient reading accuracy for the results to be considered reliable. Based on the performed experiments, the expressions for creep modulus, relaxation modulus and temperature shift factor were determined. After such analytical expressions in the mathematical form of the Prony series with knowledge of other physical characteristics of the material, numerical simulations were performed.

Determination of Creep Modulus from Creep Test
The well-known Heaviside function shown in Figure 6a is used to determine the creep modulus in this example. As a widely used function for a mechanical test, there are no particular difficulties in controlling the input, in this case the stress (force) since it is a creep test. The result is a deformation (creep function) during the test time ε(t) (Figure 6b). tics of the material, numerical simulations were performed.

Determination of Creep Modulus from Creep Test
The well-known Heaviside function shown in Figure 6a is used to determine the creep modulus in this example. As a widely used function for a mechanical test, there are no particular difficulties in controlling the input, in this case the stress (force) since it is a creep test. The result is a deformation (creep function) during the test time ε(t) (Figure 6b).  For these tests, we also chose stress values so that the total deformation does not exceed 2% because a pronounced creep characteristic is expected in this area. In order to demonstrate the usefulness and applicability of the proposed procedure to different polymeric printed materials, we use real experimental data measured at constant load (σ o = const) for the creep test. To determine the relaxation function in this case, which is the creep modulus, Equation (5) is used, i.e., the unknown coefficients J i and τ i are determined. The nonlinear regression procedure is implemented in MATLAB software and applied here.
Material samples were prepared on the described experimental line [23] or 3D printer in a layer thickness of value d 1 = 0.09 mm. The dimensions of the samples used in these test cases are 50 × 12 × 4 mm. The results of testing and approximation are presented in Table 1, while in Figure 7 a graphical presentation of the results of approximation and experimental results is given. The number of exponential terms in the series was systematically increased from n = 1 for the roughest approximation to n = 7 for the finest. From Table 1 it can be seen that the optimization was performed when choosing the relaxation times, which means that the relaxation times are assumed to be known in order to reduce the number of unknown parameters and speed up the procedure for solving systems of nonlinear equations.
3.2 × 10 −9 3.1 × 10 −9 2.9 × 10 −9   Mechanical tests with a linear increase in load (stress/strain) over time are often used to determine creep modulus (Figure 8). This is a very demanding test and is reflected primarily in load control over time. However, this problem can be overcome quite well with the advent of the newer generation of rheometers.

Number of Terms
For this test, we selected the same 3D printed polymeric material (PLA) with samples printed in three-layer thicknesses. The aim of this test is to determine the creep modulus or more precisely to determine the coefficients in the Prony series for the creep modulus during the test time. For this test, the number of exponential terms in the Prony series was taken from N = 1, 2, 5, and 10 with the relaxation times based on previous experiments and the linear distributions over several decades. The assumption of known relaxation times is based on previous similar research, and the approach used in this study is more convenient due to the reduction in the number of unknowns and equations in the system of nonlinear equations, as a consequence of the application of nonlinear regression. The coefficients in Table 1 were applied in the expression for the creep modulus and compared with the experimentally determined values. Figure 7 for the PLA material and layer thickness d = 0.09 mm shows different deviations of the fitted parameter values from the experimental ones for different numbers of exponential members in the creep module. In the initial phase of the test, a larger deviation of the fitted ones from the experimental values was noticed for the number of members n = 1 and n = 2, while for n = 6 and n = 7 a larger deviation was noted at the end of the test. The best results are obtained for the number of exponential terms n = 3 where the maximum deviation of experimental and fitted results is up to 1.85%. The maximum deviation for n = 1 at the beginning of testing is 2.9% while for n = 7 it is 5.3%.
In Table 2 the values of the average deviation between the fitted and experimental data for the whole test period are shown. The results show that the best solution is the Prony series with n = 7 exponential terms. Mechanical tests with a linear increase in load (stress/strain) over time are often used to determine creep modulus (Figure 8). This is a very demanding test and is reflected primarily in load control over time. However, this problem can be overcome quite well with the advent of the newer generation of rheometers.
For this test, we selected the same 3D printed polymeric material (PLA) with samples printed in three-layer thicknesses. The aim of this test is to determine the creep modulus or more precisely to determine the coefficients in the Prony series for the creep modulus during the test time. For this test, the number of exponential terms in the Prony series was taken from n = 1, 2, 5, and 10 with the relaxation times based on previous experiments and the linear distributions over several decades. The assumption of known relaxation times is based on previous similar research, and the approach used in this study is more convenient due to the reduction in the number of unknowns and equations in the system of nonlinear equations, as a consequence of the application of nonlinear regression.   Figure 8b shows the experimental results for deformation over time as a consequence of the applied load for one layer (d = 0.19 mm) of the 3D printed samples. All results in the form of certain coefficients were obtained using the appropriate Equation (5) in relation to the experimental results.
These coefficients are shown in Table 3. Figure 9 shows the good agreement of the results layer with regard to the d = 0.19 mm layer of 3D printed ABS. It is also noticeable that the results match throughout the test time, unlike the previous classic creep test with a constant stress value. In these tests, the logarithmic distribution of relaxation times in individual terms of the Prony series was used, which gave quite good fitting results even for a small number of exponential terms. Unlike the previous test, a certain consistency is shown here in terms of reducing the mean deviation of the fitted from the experimental results with an increase in the number of exponential terms. All these deviations are less than 2%, which can be assessed as quite acceptable.    Figure 8b shows the experimental results for deformation over time as a consequence of the applied load for one layer (d = 0.19 mm) of the 3D printed samples. All results in the form of certain coefficients were obtained using the appropriate Equation (5) in relation to the experimental results.
These coefficients are shown in Table 3. Figure 9 shows the good agreement of the results layer with regard to the d = 0.19 mm layer of 3D printed ABS. It is also noticeable that the results match throughout the test time, unlike the previous classic creep test with a constant stress value. In these tests, the logarithmic distribution of relaxation times in individual terms of the Prony series was used, which gave quite good fitting results even for a small number of exponential terms. Unlike the previous test, a certain consistency is shown here in terms of reducing the mean deviation of the fitted from the experimental results with an increase in the number of exponential terms. All these deviations are less than 2%, which can be assessed as quite acceptable.

Determination of Relaxation Modulus from the Constant Strain-Rate Test
For this test, we selected identical 3D printed polymeric material (ABS) from the previous Section 3.1.2 with samples printed in layer thickness d = 0.19. The aim of this test is to determine the relaxation modulus or more precisely to determine the coefficients in the Prony series for the relaxation modulus (E(t)) during the test time.
As in the previous test, the number of exponential members in the Prony series was taken from N = 1, 2, 5, and 10 ( Table 4), but this time the relaxation times were taken from previous research. Figure 10 shows the method of loading 3D printed samples of polymeric material and shows the linear stress ratio over time. The dimensions and geometry of the samples are identical to those in Sections 3.1.1. and 3.1.2.

Determination of Relaxation Modulus from the Constant Strain-Rate Test
For this test, we selected identical 3D printed polymeric material (ABS) from the previous Section 3.1.2 with samples printed in layer thickness d = 0.19. The aim of this test is to determine the relaxation modulus or more precisely to determine the coefficients in the Prony series for the relaxation modulus (E(t)) during the test time.
As in the previous test, the number of exponential members in the Prony series was taken from n = 1, 2, 5, and 10 ( Table 4), but this time the relaxation times were taken from previous research. Figure 10 shows the method of loading 3D printed samples of polymeric material and shows the linear stress ratio over time. The dimensions and geometry of the samples are identical to those in Section 3.1.1. and Section 3.1.2. to determine the relaxation modulus or more precisely to determine the coefficients in the Prony series for the relaxation modulus (E(t)) during the test time.
As in the previous test, the number of exponential members in the Prony series was taken from N = 1, 2, 5, and 10 ( Table 4), but this time the relaxation times were taken from previous research. Figure 10 shows the method of loading 3D printed samples of polymeric material and shows the linear stress ratio over time. The dimensions and geometry of the samples are identical to those in Sections 3.1.1. and 3.1.2.    Figure 11 shows the experimental results for stress over time as a consequence of the applied load for layer thickness d = 0.19 mm of 3D printed samples. All results in the form of certain coefficients were obtained using the appropriate expression (5) for stress (σ(t)) in relation to the experimental results.
Similar to the test case of a constant increase in stress, the use of this test allowed a fairly precise determination of the relaxation function even with a small number of members in the Prony series. The previously adopted logarithmic distribution of relaxation times and the use of the derived expression (5) gave deviation values of fitted and experimental results that are less than 2%.

Determination of Relaxation Module from Dynamic Test "Frequency Sweep"
During this dynamic test, the frequency changes, while the amplitude of deformation or the amplitude of the shear stress is kept constant. Data at low frequencies describe the behavior of samples at slow load changes. Analogous to this, the fast load behavior is expressed at high frequencies. For this example, ABS material was taken from which test samples were made and tested on a rheometer. The relaxation module was presented in the form of a Prony series where the number of members gradually varied towards higher values. The results are presented in the following Figure 12 and Table 5.  Figure 11 shows the experimental results for stress over time as a consequence of the applied load for layer thickness d = 0.19 mm of 3D printed samples. All results in the form of certain coefficients were obtained using the appropriate expression (5) for stress (σ(t)) in relation to the experimental results. Similar to the test case of a constant increase in stress, the use of this test allowed a fairly precise determination of the relaxation function even with a small number of members in the Prony series. The previously adopted logarithmic distribution of relaxation times and the use of the derived expression (5) gave deviation values of fitted and experimental results that are less than 2%.
3.1.4. Determination of Relaxation Module from Dynamic Test "Frequency Sweep" To approximate the mentioned experimental data, expression (8) with the previously mentioned nonlinear regression was used. From the presented Figure 12 we have a good agreement of experimental and approximated results. The results can also be presented in the following tabular form: For this procedure, the maximum deviation of the measured and approximated value of the shear stress storage module was calculated. Figure 12b shows the decrease in the maximum deviation of the two compared sizes with the increase in the number of members in the Prony series as expected. It can also be seen that in the first approximation with the Prony series in which we have only one member that the deviation is less than 1%, which is quite acceptable.

Determination of Time-Temperature Coefficient aT
To calculate the shift factor, a series of measurements of the relaxation module at different temperatures starting from T1 = 30 °C to T4 = 60 °C with an increase of ΔT = 10 °C was performed (Figure 13). The value of the so-called experimental window EW = 5.0 × 10 +1 s. The above expression (18) was used primarily to determine the shift factor for an individual temperature difference (e.g., 30 °C and 40 °C).
To approximate the mentioned experimental data, expression (8) with the previously mentioned nonlinear regression was used. From the presented Figure 12 we have a good agreement of experimental and approximated results. The results can also be presented in the following tabular form: For this procedure, the maximum deviation of the measured and approximated value of the shear stress storage module was calculated. Figure 12b shows the decrease in the maximum deviation of the two compared sizes with the increase in the number of members in the Prony series as expected. It can also be seen that in the first approximation with the Prony series in which we have only one member that the deviation is less than 1%, which is quite acceptable.

Determination of Time-Temperature Coefficient a T
To calculate the shift factor, a series of measurements of the relaxation module at different temperatures starting from T 1 = 30 • C to T 4 = 60 • C with an increase of ∆T = 10 • C was performed ( Figure 13). The value of the so-called experimental window EW = 5.0 × 10 +1 s. The above expression (18) was used primarily to determine the shift factor for an individual temperature difference (e.g., 30 • C and 40 • C).
To determine the value of the shift factor between the two isothermal values, a nonlinear regression implemented in MATLAB software was used again, which was run by the "lsqfit" command. The procedure was further repeated for other values of the shift factor between the values of the relaxation module for the reference temperature (in this case T ref = 30 • C) and the values of the relaxation module for the experimental isothermal conditions (T 2 = 40 • C, T 3 = 50 • C and T 4 = 60 • C). After determining the individual values of the shift factor for all temperature values in relation to the reference, expression (19) and nonlinear regression were used to determine the unknown constants C 1 and C 2 . Using the mentioned methodology and experimental data, the values of WLF constants C 1 = 0.01593 and C 2 = 28.13 were obtained for the reference temperature T ref = 30 • C. The values for the time-temperature factor are presented graphically in Figure 13b. A good agreement between the experimental (measured) results and the results obtained using nonlinear regression was found. The MATLAB software package and the "lsqfit" command were used to approximate the experimental results. To determine the value of the shift factor between the two isothermal values, a nonlinear regression implemented in MATLAB software was used again, which was run by the "lsqfit" command. The procedure was further repeated for other values of the shift factor between the values of the relaxation module for the reference temperature (in this case Tref = 30 °C) and the values of the relaxation module for the experimental isothermal conditions (T2 = 40 °C, T3 = 50 °C and T4 = 60 °C). After determining the individual values of the shift factor for all temperature values in relation to the reference, expression (19) and nonlinear regression were used to determine the unknown constants C1 and C2. Using the mentioned methodology and experimental data, the values of WLF constants C1 = 0.01593 and C2 = 28.13 were obtained for the reference temperature Tref = 30 °C. The values for the time-temperature factor are presented graphically in Figure 13b. A good agreement between the experimental (measured) results and the results obtained using nonlinear regression was found. The MATLAB software package and the "lsqfit" command were used to approximate the experimental results.

Examples of Numerical Simulations
Numerical calculations were performed for polymeric materials for which the relaxation modulus is already known as well as for the material for which the relaxation modulus was determined in a previous experimental analysis. Calculations were performed for pure mechanical load, temperature influence as well as load with simultaneous cooling of the material. Experimental data or available analytical solutions were used to compare the results.

Relaxation Test
In this example, the relaxation test for the sample in Figure 14 made of a viscoelastic material is considered. The horizontal movement of the left side of the rod is limited, while the uniaxial load is set on the right side of the rod. The resolution domain is divided into 10 CV with a time step δt = 1 s.

Examples of Numerical Simulations
Numerical calculations were performed for polymeric materials for which the relaxation modulus is already known as well as for the material for which the relaxation modulus was determined in a previous experimental analysis. Calculations were performed for pure mechanical load, temperature influence as well as load with simultaneous cooling of the material. Experimental data or available analytical solutions were used to compare the results.

Relaxation Test
In this example, the relaxation test for the sample in Figure 14 made of a viscoelastic material is considered. The horizontal movement of the left side of the rod is limited, while the uniaxial load is set on the right side of the rod. The resolution domain is divided into 10 CV with a time step δt = 1 s. The relaxation modulus of the material in the form of the Prony series was previously experimentally determined and is presented in Table 3. The loads were performed on the right side of the rod by specifying the displacement increment ( Figure 14) which can be described in mathematical form as follows: In the given example, 2D spatial discretization with a uniform orthogonal numerical mesh was used, although the example can also be solved in the 1D domain. For regions 2 and 4 the displacement vector is known at the boundary (Equation (22)), for region 1 (Fig-Figure 14. The relaxation modulus of the material in the form of the Prony series was previously experimentally determined and is presented in Table 3. The loads were performed on the right side of the rod by specifying the displacement increment ( Figure 14) which can be described in mathematical form as follows: In the given example, 2D spatial discretization with a uniform orthogonal numerical mesh was used, although the example can also be solved in the 1D domain. For regions 2 and 4 the displacement vector is known at the boundary (Equation (22)), for region 1 ( Figure 15) the boundary condition is the symmetry plan (Equation (23)), while for region 3 the known (free surface) load is valid, i.e., the boundary condition is also from Equation (22). The relaxation modulus of the material in the form of the Prony series was previously experimentally determined and is presented in Table 3. The loads were performed on the right side of the rod by specifying the displacement increment ( Figure 14) which can be described in mathematical form as follows: ε = 5.75 • 10 ; t ≥ 0 . (30) In the given example, 2D spatial discretization with a uniform orthogonal numerical mesh was used, although the example can also be solved in the 1D domain. For regions 2 and 4 the displacement vector is known at the boundary (Equation (22)), for region 1 (Figure 15) the boundary condition is the symmetry plan (Equation (23)), while for region 3 the known (free surface) load is valid, i.e., the boundary condition is also from Equation (22).
To solve this case, a grid with (50 × 20) 1000 control volumes was used. Figure 16 shows that as the number of exponential members in the relaxation module increases, the numerical results increasingly agree with the experimental results.   To solve this case, a grid with (50 × 20) 1000 control volumes was used. Figure 16 shows that as the number of exponential members in the relaxation module increases, the numerical results increasingly agree with the experimental results.

Creep Recovery Test
For the next test, the same material (ABS) was selected from the relaxation test (3.2.1) with identical characteristics of the relaxation modulus, while the load consists of one period of constant load (σ o = 7.5 MPa) to the moment t = 100 s, followed by unloading (σ o = 0 MPa) ( Figure 17). The geometry and domain of the solution are identical to those in Section 3.2.1, and the load limit condition was applied to region 2 of Figure 15. Other boundary conditions as well as geometry are the same as in Figure 14a. An expression with one free and one exponential term according to Table 4 was used for the relaxation module. The relaxation modulus of the material in the form of the Prony series was previously experimentally determined and is presented in Table 3. The loads were performed on the right side of the rod by specifying the displacement increment ( Figure 14) which can be described in mathematical form as follows: In the given example, 2D spatial discretization with a uniform orthogonal numerical mesh was used, although the example can also be solved in the 1D domain. For regions 2 and 4 the displacement vector is known at the boundary (Equation (22)), for region 1 (Figure 15) the boundary condition is the symmetry plan (Equation (23)), while for region 3 the known (free surface) load is valid, i.e., the boundary condition is also from Equation (22).
To solve this case, a grid with (50 × 20) 1000 control volumes was used. Figure 16 shows that as the number of exponential members in the relaxation module increases, the numerical results increasingly agree with the experimental results.   For this case, an analytical solution can be found that takes the form: where the coefficients in the creep modulus are determined on the basis of a three-parameter model: where r is the relaxation time in the creep modulus. More about performing an analytical solution can be found in the literature [11]. For the next test, the same material (ABS) was selected from the relaxation test (3.2.1) with identical characteristics of the relaxation modulus, while the load consists of one period of constant load (σo = 7.5 MPa) to the moment t = 100 s, followed by unloading (σo = 0 MPa) ( Figure 17). The geometry and domain of the solution are identical to those in Section 3.2.1, and the load limit condition was applied to region 2 of Figure 15. Other boundary conditions as well as geometry are the same as in Figure 14a. An expression with one free and one exponential term according to Table 4 was used for the relaxation module. For this case, an analytical solution can be found that takes the form: where the coefficients in the creep modulus are determined on the basis of a three-parameter model: where r is the relaxation time in the creep modulus. More about performing an analytical solution can be found in the literature [11]. In Figure 17b a fairly good agreement of the numerical results with the measured experimental results (−o) is shown. In the same figure it can be seen that the creep phenomenon is most pronounced at the beginning of the test.

Plane-Strain Time-Varying Test
In this case, the material with the experimentally determined relaxation modulus from Table 3 was also used. The relaxation module is characterized as a Prony series with two exponential terms and one free term. The sample density is known and is 900 kg/m 3 . Poisson's ratio for ABS is taken from the literature (υ = 0.35). The geometry and boundary conditions are given in the following Figure 18a. This test case is analyzed in the 2D domain, and since no analytical solution is known, the obtained simulation results are compared with the results of a numerical simulation carried out using the finite element method. The load at the boundary of the solution domain consists of a time-varying strain, while the total test time is divided into four regions. In Figure 17b a fairly good agreement of the numerical results with the measured experimental results (−o) is shown. In the same figure it can be seen that the creep phenomenon is most pronounced at the beginning of the test.

Plane-Strain Time-Varying Test
In this case, the material with the experimentally determined relaxation modulus from Table 3 was also used. The relaxation module is characterized as a Prony series with two exponential terms and one free term. The sample density is known and is 900 kg/m 3 . Poisson's ratio for ABS is taken from the literature (υ = 0.35). The geometry and boundary conditions are given in the following Figure 18a. This test case is analyzed in the 2D domain, and since no analytical solution is known, the obtained simulation results are compared with the results of a numerical simulation carried out using the finite element method. The load at the boundary of the solution domain consists of a time-varying strain, while the total test time is divided into four regions. The first region consists of a linear increase in deformation to a value of 1% (εo = 0.01, t1 = 5 s), the second is the maintenance of constant deformation until t2 = 50 s, in the third region a linear unloading is performed until t3 = 55 s (εo = 0), and in the last region the relaxation effect is monitored until t4 = 100 s. The time course of the load is also presented in Figure 18b.
An orthogonal numerical grid of 100 CVs was used to discretize the solution domain (Figure 18a), while a step of t = 1 s was chosen for time discretization. Figure 19 presents relaxation effect is monitored until t 4 = 100 s. The time course of the load is also presented in Figure 18b.
An orthogonal numerical grid of 100 CVs was used to discretize the solution domain (Figure 18a), while a step of t = 1 s was chosen for time discretization. Figure 19 presents the numerical results obtained using two numerical methods (Finite Volume Method and Finite Element Analysis). Two numerical methods were used due to the very difficult performance of an experiment of this type, and an analytical solution for this test case (υ = 0) is not available. The application of the numerical finite element method was performed using the well-known ABAQUS software with built-in subroutines for calculating viscoelastic deformations and stresses. The numerical solutions agree very well in all load regions and good agreement is especially visible at linear loading or unloading. The first region consists of a linear increase in deformation to a value of 1% (εo = 0.01, t1 = 5 s), the second is the maintenance of constant deformation until t2 = 50 s, in the third region a linear unloading is performed until t3 = 55 s (εo = 0), and in the last region the relaxation effect is monitored until t4 = 100 s. The time course of the load is also presented in Figure 18b.
An orthogonal numerical grid of 100 CVs was used to discretize the solution domain (Figure 18a), while a step of t = 1 s was chosen for time discretization. Figure 19 presents the numerical results obtained using two numerical methods (Finite Volume Method and Finite Element Analysis). Two numerical methods were used due to the very difficult performance of an experiment of this type, and an analytical solution for this test case (υ ≠ 0) is not available. The application of the numerical finite element method was performed using the well-known ABAQUS software with built-in subroutines for calculating viscoelastic deformations and stresses. The numerical solutions agree very well in all load regions and good agreement is especially visible at linear loading or unloading.

Analysis of Viscoelastic Material with Simultaneous Mechanical and Thermal Loading
For this case, the material that was previously subjected to experimental analysis in determining the relaxation modulus in Section 3.1.3 was analyzed. This is a very convenient test example because the numerically obtained data can be further directly compared with the experimental data. The thermal effect is taken into account through the use of shift factors in the form of a WLF function. The values for the coefficients in the WLF function were taken from Section 3.1.4 since the same material was used. Other physical characteristics are presented in Table 6, below.

Analysis of Viscoelastic Material with Simultaneous Mechanical and Thermal Loading
For this case, the material that was previously subjected to experimental analysis in determining the relaxation modulus in Section 3.1.3 was analyzed. This is a very convenient test example because the numerically obtained data can be further directly compared with the experimental data. The thermal effect is taken into account through the use of shift factors in the form of a WLF function. The values for the coefficients in the WLF function were taken from Section 3.1.4 since the same material was used. Other physical characteristics are presented in Table 6, below. Geometrical characteristics and boundary conditions are presented in Figure 20 where the initial temperature T 1 = 30 • C at all points of the domain.
lation to the experimentally determined results from moving the end of the rod.    The material is gradually heated by means of a linearly changed temperature T o . The design of the test device is such that a sample is fixed at one end in the experimental chamber while a constant stress load is set at the other end. The temperature rate is constant and after t = 180 s it is 60 • C. The load at the end of the rod is σ = 1.51× 10 +5 Pa.
After the geometry was made and the initial and boundary conditions were applied, the following results were obtained as in Figure 21. The presentation of the results shows a fairly good agreement between the results obtained using the numerical method in relation to the experimentally determined results from moving the end of the rod.
Geometrical characteristics and boundary conditions are presented in Figure 20 where the initial temperature T1 = 30 °C at all points of the domain.
The material is gradually heated by means of a linearly changed temperature To. The design of the test device is such that a sample is fixed at one end in the experimental chamber while a constant stress load is set at the other end. The temperature rate is constant and after t = 180 s it is 60 °C. The load at the end of the rod is σ = 1.51× 10 +5 Pa.
After the geometry was made and the initial and boundary conditions were applied, the following results were obtained as in Figure 21. The presentation of the results shows a fairly good agreement between the results obtained using the numerical method in relation to the experimentally determined results from moving the end of the rod.

Discussion
The experimental research of relaxation functions presented in this paper can be applied to a wide range of polymeric materials, starting from "pure" plastics, such as ABS, to materials formed as a mixture of polymers and some basic material (e.g., wood or a combination of wood and glue). The selected application examples demonstrate the specific determination of the coefficients in the Prony series based on the experimental values of the relaxation function. Two examples of determining the analytical form of the relaxation function show that the methodology for determining the unknown coefficients in the Prony series is applicable in the real time (t) domain as well as in the angular velocity (ω)-domain. After determining the analytical expression under isothermal conditions, the influence of temperature on the values of the relaxation function was investigated, where the time-temperature principle was applied. The application of this principle was performed within the experimental window, where the "shift-factor" values were determined using the WLF equation and nonlinear regression between the experimental values and the assumed analytical expression. After this analysis, a numerical calculation was performed for non-isothermal conditions with such an experimentally determined relaxation function and shift factor, with a small difference between the numerical and experimental results. A rheometer test device was used for testing the samples, as well as a universal ZwickRoell compression and extension testing machine. Both devices provided quite good control of the input parameters, primarily the control of the sample load, whether it was a creep or relaxation test. Nonlinear regression, i.e., the well-known least squares method using the subroutines in MATLAB, very quickly gave good results of "fitting" the experimental data for an arbitrary and optimized set of coefficients in the Prony series. For the application of the recursive algorithm, the finite volume method was used, which was implemented in earlier papers for the problems of viscoelastic solids. The basis for all numerical calculations is the software package SA-2D developed for the analysis of thermo-viscoelastic deformations. Finally, it was shown that for polymer 3D printed materials a relaxation function can be determined in the form of a Prony series based on the experiments performed in this study. Based on the analytical expressions of the Prony series and shift factors, the numerical calculations yielded results that differed very little from the experimental ones. The numerical results using the recursive algorithm did not deviate significantly from the experimental ones even when considering problems with non-isothermal conditions.