Modeling of thermo‐viscoelastic material behavior of glass over a wide temperature range in glass compression molding

In glass compression molding, most current modeling approaches of temperature‐dependent viscoelastic behavior of glass materials are restricted to thermo‐rheologically simple assumption. This research conducts a detailed study and demonstrates that this assumption, however, is not adequate for glass molding simulations over a wide range of molding temperatures. In this paper, we introduce a new method that eliminates the prerequisite of relaxation functions and shift factors for modeling of the thermo‐viscoelastic material behavior. More specifically, the temperature effect is directly incorporated into each parameter of the mechanical model. The mechanical model parameters are derived from creep displacements using uniaxial compression experiments. Validations of the proposed method are conducted for three different glass categories, including borosilicate, aluminosilicate, and chalcogenide glasses. Excellent agreement between the creep experiments and simulation results is found in all glasses over long pressing time up to 900 seconds and a large temperature range that corresponds to the glass viscosity of log (η) = 9.5 – 6.8 Pas. The method eventually promises an enhancement of the glass molding simulation.

down to a solid state. Depending on the temperature difference between the glass and mold during pressing, the glass molding process is distinguished by isothermal glass molding or precision glass molding (PGM) and nonisothermal glass molding (NGM). 2 Regardless, the glass undergoes the heating stage to a liquid-like state at nearly softening temperature (S p ) in PGM or even above S p in NGM and the following pressing and cooling stages to a solid-like state, typically below transition temperature (T g ), where its desirable shape is mostly achieved.
In this temperature regime, viscoelastic properties dominate the deformation of the glass under thermo-mechanical loads. Previous studies demonstrate essential contributions of viscoelastic properties to the final shape of molded lenses. [3][4][5] Numerical predictions of the final lens shape requires a coupling thermo-mechanical constitutive model to comprehend thermo-viscoelastic responses of glass deformation over the entire molding process. For mechanical modeling of the viscoelastic responses, a large number of rheological models have been investigated in the field of glass molding simulation over years. In the early days, behavior of glass at the molding temperature was commonly considered the pure viscous flow, where flow stress is proportional to the strain rate. [6][7][8] Other studies use either simple Maxwell or Kelvin model to approximate the stress relaxation and creep responses, respectively. 9,10 Those simple models are chosen to qualitatively study the form deviation of molded glass shapes and stress characteristics. A later study by Zhou et al reports that neither Maxwell nor Kelvin model satisfactorily describes the experiment data in the temperature regime of glass molding, but Burgers model-a combination of both Maxwell and Kelvin model-does. 11 The Burgers model, in fact, can be represented by a twoterm generalized Maxwell model with an identical mathematical description. 12 The generalized Maxwell model, a combination of several Maxwell elements in parallel, has been widely utilized in many recent studies to characterize the viscoelastic property of glass above the transition temperature. [13][14][15] The sufficient number of Maxwell elements is determined by fitting the relaxation data, either gained by stress relaxation experiments or derived from creep experiments. It is emphasized that the number of Maxwell elements strongly depends on the temperature chosen for the experimental characterization. For example, a six-term generalized Maxwell model is commonly found for investigating temperatures in the vicinity of T g , 16 while three terms are used at higher characterizing temperatures. 17 In addition, the thermal-mechanical modeling of the viscoelastic responses necessitates an incorporation of the temperature effect into the mechanical models. Most of the present studies relied on thermo-rheologically simple (TRS) modeling approach in the field of glass molding simulation. [18][19][20] The TRS assumption is widely employed due to the ease of numerical implementation where the coupling temperature-mechanical model can be simplified as a solely mechanical model by incorporating the so-called reduced time, defined by a shift function. It is emphasized that this implementation is valid if and only if the shape of stress relaxation functions at every temperature, when plotted against a logarithmic time, is fully identical. However, Scherer stresses that, albeit a good approximation near the glass transition, it fails at high temperature where the glass exhibits liquid-like viscoelasticity. 21 In particular, sensitivity studies of the TRS assumption highlight that non-negligible deviation occurs in the numerical predictions of the molded lens shape 22 or creep displacements. 5,23 Finally, to characterize the thermo-viscoelastic material behavior of the glass, the most common approach, saying conventional approach, is to perform experiments, either creep or stress relaxation, at several temperatures. The stress relaxation experiment is likely more favorable because the model's parameters, for example, the generalized Maxwell model, can be directly obtained by fitting the stress relaxation data measured at one specific temperature, called the reference temperature. 24 The disadvantage of this approach is that the stress relaxation experiments are mostly conducted at temperatures in the vicinity of T g , which are quite far from the actual molding temperature, normally near S p . In this context, the viscoelastic properties at the actual molding temperature are determined by using the shift function. In contrast, the characterization of viscoelastic properties at the real molding temperature is restricted to the creep experiment. At such high temperature, in fact, the glass becomes liquid-like where its relaxation time is extremely fast (τ«seconds) so that the data obtained from the stress relaxation experiment is not sufficient and reliable for postprocessing. Based on the creep approach, the mold displacements are acquired and used to compute creep compliances. 25,26 The drawback of this approach is that the temperature-dependent viscoelasticity is not accomplished. Instead, the modulus relaxations need to be derived from the creep compliance functions using interconversion theory. Relying on the creep setup, in addition, previous studies realize a non-negligible discrepancy between the experiment and numerical prediction of the creep displacements. 26,27 According to their thoughts, the discrepancy is attributed to friction at the glass-mold contact interface. Nevertheless, satisfactory verifications of the creep displacement cannot be obtained by varying the friction coefficients. A recent study by Yu et al shows that if the frictional disturbance can be diminished from the creep data so that mere purity of the viscoelastic property is secured by creep experiments, good agreement can be found. 28 Three major contributions are aimed in this research. First, a comprehensive study on the conventional characterization of the thermo-viscoelastic properties derived from the creep experiments is presented. The study aims to demonstrate the rationale for the mismatches between the simulation and creep | 2793 experiments reported in earlier studies; and the following is to provide a solution to such existing problem. Furthermore, the thermo-mechanical modeling using TRS over a wide range of molding temperature is examined. The inadequacy of TRS over the wide temperature range is then elaborated. Finally, this research introduces a new thermo-viscoelastic modeling approach to overcome this deficit. According to the proposed method, the prerequisite of the relaxation functions and shift factors by using the TRS assumption to characterize the thermo-viscoelastic properties can be completely eliminated. Experimental validations are carried out for three different glass types commonly used in compression molding for glass optic production. Excellent agreement between experiment and simulation results found in all three glasses strengthens the validity of our method.

THERMO-VISCOELASTICIT Y
The constitutive thermo-rheological model of a viscoelastic material determines the relationship among stress, strain, strain rate, and temperature. In view of mechanical modeling, the stress-strain relation can be either described by differential forms or hereditary integrals. 29 The integral representation is commonly chosen as its convenience for incorporating temperature into the mechanical model and flexibility to manipulate the measured viscoelastic properties, ie, creep and stress relaxation. For linear viscoelastic theory, the relation between stress and strain is linear and time dependent, represented by the following three-dimensional (3D) integral forms, using Boltzmann superposition principle: where σ and ε are stress and strain tensor, respectively; and are relaxation modulus-and creep compliance tensors, respectively; t is the actual observation time, and t′ is an arbitrary time prior to t.
The rate of relaxation rapidly varies with the change of temperature. For incorporating the temperature into the mechanical formulation, the conventional approach relies on the concept of time-temperature superposition principle, mainly based on the TRS property. Based on this concept, a single parameter is introduced to combine the effects due to both time and temperature. Then, Equations (1) and (2) in varying temperature becomes: where is the so-called "reduced time", which is related to the real time t by a temperature shift factor a(T) and a reference temperature T r . In case of transient temperature conditions, ie, temperature varying with time, the reduced time is defined as follows: where the shift factor is given by: In Equation (6), τ is the relaxation time. Here a(T) defines the horizontal shift from the relaxation curve of the computed temperature Tto that of the reference temperature T r , called the master curve.
As mentioned earlier, characterization of the viscoelastic responses can be performed by either the creep or stress relaxation experiments. In fact, both these phenomena are the viscoelastic properties of the same material; therefore, they can be interconverted. According to the interconversion theory, the relation between creep compliance and relaxation modulus can be defined as a compact form in s-domain by applying Laplace transform: Such a relation can be written in the time domain by performing the inverse Laplace transform: Using Equation (7), the relaxation modulus is derived from the creep compliance on the s-domain, and following the relaxation modulus on the time domain (t) is defined by the inverse Laplace transform, using Equation (8). Then, considering the TRS assumption, the shift factor is determined from the modulus relaxations similar as the stress relaxation approach.
Finally, in 3D formulation of the viscoelastic responses, the stress--strain relation, eg, Equation (1), is commonly expressed in a form where the deviatoric and volumetric parts are separated. In this regard, the shear and bulk (volumetric) relaxation moduli are employed instead of Young modulus. Instantaneous shear and bulk moduli can be computed from Young modulus and Poisson ratio. In the viscoelastic regime, shear deformation is typically more dominant than bulk deformation 30 ; hence, time-independent bulk modulus, ie, no relaxation, is commonly agreed. Accordingly, the shear relaxation modulus and its temperature-dependent characteristic are the most essential requirements for the characterization of the thermo-viscoelastic properties.

| Experimental procedure
The characterization of the glass viscoelastic properties is based on the creep approach using the uniaxial compression setup. The creep experiments are conducted in a precision glass-molding machine Toshiba GMP 207HV. Borosilicate SUPRAX ® 8488, a "long" glass type, is chosen for the investigation of the thermo-viscoelastic responses over a wide temperature range. Thermal and mechanical properties of this glass are provided in Table 1. The glass specimen has a cylinder shape. In addition, a high-temperature resistance glassy-carbon SIGRADUR ® is used as the mold material for two essential purposes. First, this mold material allows the molding experiments at high temperatures, even above S p , without any glass sticking problem. Second, by fine polishing, the surfaces of the mold pair is mirror-like with the roughness of R a = 2 nm. 31 Such ultraprecision mold surface enables a small friction coefficient at the glass-mold interface, below 0.1. 32 As indicated in the previous works, 5,26,28 the low friction coefficient is essential to minimize the error that leads to the impurity of the viscoelastic properties derived from the creep experiments due to the frictional deformation of the glass. Figure 1 exhibits the polished glassy-carbon molds and the glass specimen used for the creep experiments. Besides, demonstrations of pressed glasses at high molding temperatures in the vicinity of S p are presented. Creep experiments are performed at four specific temperatures ranging from 650°C to 790°C, which corresponds to the viscosity range η = 10 9.5 -10 6.8 Pas. Such viscosity range, from slightly above T g to nearly S p , is of interest for this study because this is a typical molding temperature range for PGM. Besides, the glass deformation from a preform to nearly final shape is mostly occurred in this temperature range by NGM, 33,34 before it is sent to a separate annealing oven. Please note that the experimental characterization of the viscoelastic properties is conducted for equilibrium glasses. The glass specimens are slowly heated up with a rate of 5 K/min and then undergone a sufficient long soaking, approximately 20 minutes, in order to erase the thermal history of the specimens before the creep experiments. Within the chosen temperature range, this soaking time is expected to satisfactorily equilibrate the glass specimens. For the creep compression, a pressing time of 900 seconds is held for all experiments. The mold displacement is recorded and used in the following steps to derive the thermo-viscoelastic properties of the glass. Figure 2A introduces the mold displacements of four testing temperatures.

| Computational procedure of creep compliance
The computational procedure to obtain the viscoelastic properties from the creep displacements has been presented in the earlier works. 25,26 For the completeness of the paper, most important aspects are highlighted. First, true stress and true strain are taken to account for the large deformation of the glass specimens. By using the recording displacement history u(t) and the initial length of the specimen L 0 , the true strain ε(t) is given as follows:

T A B L E 1 Material properties of glass
Following the computation of the true strain, the true stress σ(t) is calculated by using the initial cross-section A 0 of the specimen and applying force F: Here, please note that the applying force F = 250N is remained for all experiments of the entire temperature range. This deployment is aimed to avoid the effect due to the pressure-dependent viscosity, 35 which is not considered in this study.
The true stress is calculated by considering a constant volume of the deformed glass specimen. Theoretically, the creep test requires an instantaneous applying stress σ 0 that holds constantly during the entire experiment. Under experimental conditions, however, such requirement cannot be fulfilled. Since the machine solely controls the applying force, the significant increase of the cross-section during the deformation of the glass specimen results in a decrease of the compression stress. To account for the change of the true stress as a variable with time, σ(t), the stress input is approximated by the sum of serial constant stress inputs σ i . 29 According to the Boltzmann superposition principle, the strain output under the variable stress input equals to the sum of all discretized strain outputs i resulting from each corresponding discretized stress inputs σ i . Now, it is important to point out that only vertical displacements of the mold, ie, one-dimensional (1D) data are acquired by the uniaxial compression setup. Thus, Equation (2) needs to be rewritten in the 1D form. In addition, by separating the initial applied stress (σ(0) = σ 0 ) and the variable stress, for the sake of creep compliance computation, an alternative form of Equation (2) becomes: This convolution integral can be solved by a finite difference scheme using discretized representations of the true stress σ n , true strain ε i , and creep compliance J i , defined as follows: Rearranging Equation (12), the discretized creep compliance J i is given as follows: Thus, from the recording displacement inputs and the initial geometry of the glass specimens, using Equations (9-10), the creep compliances J(t) are achieved by solving Equation (13) iteratively. Figure 2B provides the creep compliances computed from the displacement inputs u(t).

| Determination of relaxation modulus
The relaxation modulus can be achieved if the creep compliance function is known, according to Equation (7). Since the creep compliance contains discrete data, the creep function is realized by curve fitting. The fitting of the creep compliance data requires a material model. For example, generalized Kelvin with three terms is used in previous studies. 17,26 For the liquid-like creep behavior characterized at high temperature range investigated in this study, instead, the 4 element-Burgers model is employed. The phenomenological representation of the Burgers model is shown in Figure 3A. The creep function J(t), accordingly, is defined as follows: where E 1 , E 2 , η 1 and η 2 are the fitting parameters of the Burgers model. There are several advantages of utilizing the Burgers model. First, as demonstrated in Figure 3B, the Burgers model illustrates an outstanding fitting accuracy over the entire temperature range of investigation. An earlier study of the author has indicated that the Burgers model precisely predicts the creep displacements of B270i glass in a wide molding temperature range, which correspond to viscosity η = 10 11.7 − 10 7 Pas. 36 Second, it is emphasized that because the Burgers model consists of a Maxwell and a Kelvin element, it is capable of representing both creep and stress relaxation phenomena. As a result, the relaxation function can be directly computed from the fitting parameters obtained from the creep compliance. This property helps to avoid the second fitting of the relaxation modulus to attain the parameters of the relaxation modulus function as pursued in most of the previous works. 25,26 Moreover, the Burgers model is mathematically equivalent to a two-term Prony series; therefore, the parameters obtained from the creep compliance function (Equation 14) can be used to compute the Prony series parameters, ie, weight factor g i and relaxation time τ i (i = 1,2), in a straightforward manner. Details of the mathematical equivalence between the Burgers model and two-term Prony series are described in Appendix A. When the Prony series parameters are known, the material model can be simply implemented into most of the commercial FEM simulation software, eg, ABAQUS, using the hereditary integral formulation. Finally, compared to generalized Maxwell model with 3-6 terms commonly employed in other works, the Burgers model contains the fewest fitting parameters (4 parameters) that bring essential advantages for the proposed modeling method of the thermo-viscoelastic material behavior presented in Chapter 4.
Since four parameters of the model are already realized by fitting the creep compliance using Equation (14), the shear parameters of the Burgers model, G 1 , G 2 , μ 1 and μ 2 , can be calculated as follows: Using the shear parameters in Equation (15A), Prony series coefficients, including weight factors g i and τ i (i = 1, 2), are achieved as demonstrated in Appendix A (Equations A6 and A9). Table 2A provides the computed coefficients of four investigating temperatures. Please note that as these coefficients are calculated directly from the creep compliance parameters, the second fitting of the relaxation modulus curves is no longer necessary. Using these coefficients, Figure 4 plots the normalized shear relaxation moduli (ψ(t) = G(t)/G 0 ) of four testing temperatures, accordingly.

| Numerical verification of shear relaxation parameters
A Finite Element (FE) model is established to verify the shear relaxation parameters derived from the uniaxial creep experiments. Process inputs and boundary conditions of the FE model are described in Figure 5. For the numerical verification of the creep displacements, ie, the glass thickness changes, only the glass specimen is considered. The glass has a cylinder shape and is modeled as an axisymmetric deformable body. Vertical displacements of bottom surface are constrained, while the top surface is coupled with the applying force F taken from the experiments. Since mold parts are not necessarily included in the FE model, friction between the glass and molds is neglected. In addition, avoidance of the dimensionality effect on the numerical verifications is particularly taken. Before each creep experiment, the initial length and diameter of the glass specimens are measured with an accuracy of 3 digits. Then, these geometrical inputs are set for the respective simulation of the same experiment temperature. The numerical verification is performed using ABAQUS/CAE 2017. The glass properties and viscoelastic parameters are given in Tables 1 and 2A. Figure 6A presents the simulation results of the creep displacements. When compared with the experiment data, however, the simulation results show obvious disagreements found in all testing temperatures. In fact, the numerical predictions of the creep displacement are always underestimated. As mentioned earlier, observations of such discrepancy were reported in the previous works. 26,27 In their studies, the mismatch is believed due to friction at the glass and mold contact surfaces. They argued that the displacement acquired by the compression setup is not purely viscoelastic creep deformation of the glass but contains a frictional resistance at the interface. The presence of friction violates the nature of the creep data used for the derivation of the viscoelastic parameters, and it causes the discrepancy. Several studies investigated the influence of the friction on the viscoelastic constants using the creep compression setup. It is demonstrated that the accuracy of the obtained viscoelastic parameters can be enhanced if the friction between the glass and mold is sufficiently small, 23 as approached in this research, or the error due to the friction can be dismissed from the displacement data. 28 The authors agree that the friction at the glass-mold interface, which cannot be avoided by the compression setup, may violate the pure creep deformation due to viscoelastic property of an investigating material. In this study, however, another reason that crucially results in the disagreement between the experiment and simulation results is elucidated. The argumentation is that if the creep displacement is used as the input data for the derivation of viscoelastic parameters, and the numerical verification-indeed a reverse computation-must bring an identical result. Furthermore, because such input displacement comprises the deformation due to friction, the friction at the glass-mold interface should not be reconsidered in the FE simulation for the verification. 25 According to the thoughts, the simulation and experiment displacements have to be identical. So why do such large discrepancies observed in Figure 6A come from? Described in Sections 3.2 and 3.3, the shear relaxation modulus is derived from the uniaxial displacement, ie, the 1D experiment data, and is computed using the 1D viscoelastic formula (Equation 15A). Viscoelastic theory performed in ABAQUS simulation platform is developed for generalization of the integral representation to three dimensions. 37 In 3D formulation, the stress tensor is split into volumetric and deviatoric parts for the ease of modeling the nature of volume-preserving materials such as viscoelasticity. Here, it is emphasized that while the volumetric (dilatational) tensor carries out the bulk or volumetric deformation of the materials, the deviatoric tensor solely accounts for the shear deformation, meaning that no volumetric changes are taken. Meanwhile, as discussed earlier, the bulk deformation is time-independent, ie, no viscoelastic property. For this reason, the Poisson's ratio given in Equation (15A) must be set to 0.5 instead. This implementation is essential to correct the computation of the shear parameters of the constitutive model given in the 3D viscoelastic formulation while taking no account of the deformation due to viscoelastic bulk. 38 The time-independent bulk modulus, in fact, results in an elastic response, 5 and therefore it can be accounted by the elastic component of the Burgers model (E 1 ). Eventually, the shear relaxation parameters are defined as shown in the following revised equation: Table 2B provides the revised Prony series coefficients calculated from shear relaxation parameters using Equation (15B). Based on these coefficients, simulations are rerun using the same setting as described in Figure 5. Figure 6B plots the simulation results and the experiment creep displacements. It is clear that great agreements of all experiments are found. Please remind the underestimation of the simulated creep displacements observed in Figure 6A. This underestimation results from the contribution of the bulk deformation implicit in the shear parameters by Equation (15A), which is certainly not when the deformation of the viscoelastic bulk is excluded from the computation of the shear parameters using Equation (15B). The validated results emphasize the essential implementation of the shear relaxation coefficients in the 3D formulation for the simulation of the glass molding processes.
Up to now, the viscoelastic parameters derived from the creep displacement at individual testing temperature are fully validated. For the glass molding simulation, however, temperature-dependent viscoelastic parameters are required. For modeling the thermo-viscoelastic responses, most of the previous studies relied on the TRS behavior of the viscoelastic materials. The next step focuses on an examination of this thermo-viscoelastic modeling approach over the temperature range chosen in this study.

| Numerical verification of thermoviscoelastic parameters using TRS
Thermo-viscoelastic modeling using TRS necessitates a shift function. The shift function is determined from the shear relaxation modulus functions varying with temperatures ( Figure 4). First, the master curve at an arbitrary experiment temperature-the (15B) , G 2 = E 2 2 (1 + 0.5) , 1 = reference temperature T r -is selected. In this study, the shear relaxation function at 680°C is considered. Following, this master curve is shifted to the relaxation curves of other temperatures along the abscissa. Relying on the Williams-Landel-Ferry (WLF) equation, the shift function is defined as follows: 39 where C 1 and C 2 are constants obtained by fitting the shifts computed from the master curve to each individual temperature. Figure 7 exhibits the curve fitted by WLF equation and the constants. Figure 8 provides simulation results of the creep displacements using temperature-dependent viscoelasticity by using the WLF shift function. At the reference temperature (T = 680°C), the experiment and simulation results are certainly validated, exactly the same as observed in Figure 6B. However, discrepancies between the simulation and experiment data occur at other temperatures. The discrepancy becomes larger and non-negligible at temperatures further from the reference temperature. This evidence indicates the inadequate implementation of the TRS behavior for a large temperature range. Consequently, numerical predictions of the molded glass characteristics, eg, form deviation or stress, are not reliable.
The underlying cause of the observed discrepancies stem from an improper assumption using the TRS over the range of studying temperature. In principle, the application of the TRS is valid if and only if the master relaxation curve and the shifted curves are fully identical, because only a horizontal shift is acted by the TRS theory. In order to examine this requisite, Figure 9A plots the relaxation functions at four testing experiments (solid lines) and the relaxation curves shifted from the master curve (dashed lines). The plots reveal that the modulus relaxation function and the shifted curve are not fully identical. While two data have a good match within the first relaxation response, ie, fast relaxation time (τ 1 ≈10 −5 − 10 −2 seconds), an obvious distinction can be recognized when one looks a bit closer on the long-term relaxation response (τ 2 ≈10 −2 -10 3 seconds) as indicated in the zoomed plot ( Figure 9B). Similar findings can be seen in the work of Zhou et al 26 and Yu et al 28 ; however, impacts of such nonidentical relaxation response in the long term on the creep displacements have not been studied. It is obvious that the first relaxation period occurs extremely short (<0.01 second) compared to the total pressing time. In contrast, the second relaxation period, though it weighs only a small proportion of the entire relaxation spectrum, holds the entire creep experiment, starting at very early stage (0.01 second) till the end of the pressing time (900 seconds). Therefore, the second relaxation stage definitely contributes a substantial dominance, while the first relaxation barely plays a role in the overall deformation of the glass in the creep experiments. According to the definition of the TRS (Equation 6), only one relaxation time-the first relaxation time 1 -is mainly involved in the determination of the shift factor because of its dominant overlapping ratio. The second relaxation times, however, have not the same temperature dependence, leading to the dissimilarities between the master and shifted curves when using the shift factors computed from the first relaxation times. Such inadequate characteristic of the TRS makes it insufficient to represent the temperature-dependent viscoelastic property over a wide temperature range. Clearly seen in Figure  8, even a small difference of relaxation response in the long-term relaxation time can bring substantial deviations in the numerical computation of creep displacements.
Several studies attempt to enhance the model predictions by not only shifting the master curve horizontally but also vertically, called Thermo-Rheologically Complex, so that the vertical differences can be reduced. 40 However, it is stressed that the vertical shift, again, is another assumption, and it is only valid if the master and shifted curves are fully vertically identical. In the following chapter, an alternative approach is introduced to improve the modeling of temperature-dependent viscoelasticity without assuming that the thermorheological material behavior is simple or complex.

VISCOELASTIC MODELING APPROACH
This chapter introduces a new modeling approach of the thermo-viscoelastic behavior of glass materials. The proposed method aims at the characterization of the temperature-dependent viscoelastic parameters without using the relaxation moduli and shift function. To characterize the viscoelastic properties at a specific temperature, ie, isothermal rheology, the earlier verifications demonstrate that the rheological Burgers model is fully adequate to represent the creep response ( Figure 6B). For incorporating temperature into the rheological model, ie, thermorheology, the central point is to pursue a unique set of Burgers model parameters, called Burgers parameters, which are dependent on temperature, and this unique parameter set is able to represent the glass viscoelastic responses over the entire temperature range. As previously presented, the relaxation parameters, ie, relaxation time τ i and weight factors g i (i = 1, 2), can be directly calculated from four Burgers parameters (Equations A6 and A9; Appendix A). It means that if there exists a temperature-dependent function for each individual Burgers parameter, the relaxation parameters at any temperature can be determined while the prerequisite of the relaxation functions and shift factors is no longer necessary.

| Determination of Burgers parameters
Our method focuses on finding a parameter set of the Burgers model, , which is able to represent both creep and stress relaxation over the entire temperature range. Using the creep compliances at four experiment temperatures ( Figure 3B), the Burgers parameters of each temperature are obtained by a nonlinear curve fitting. Figure 10 presents the fitting values of these parameters in natural logarithmic scales varying with temperatures. Here, it is noted that the first spring of the Burgers model, E 1 , represents the instantaneous viscoelastic response, which is accounted for the Young modulus of glass materials. This parameter can be temperature dependent, measured by several available techniques, such as Brillouin light scattering 3 or impulse excitation. 41 However, since most of glass molding simulations consider a constant Young modulus, usually available in glass datasheet, the parameters E 1 is assumed to be temperature independent. Thus, this value is constant given at room temperature (Table 1). The relations of other parameters, ie, E 2 , η 1 and η 2 , with temperature, plotted in Figure 10, reveal a unique trend where the properties of these parameters are decreasing with increasing temperatures.
In the next step, such relations between each Burgers parameter and temperature demand a best mathematical description, ie, a relevant temperature-dependent model. In this regard, it is necessary to figure out a realistic model using the fewest possible number of fitting parameters. Since the observing trends shown in Figure 10 are similar to that of the glass viscosity, existing models for the equilibrium viscosity of glass can be a good reference. In the transition region, most of the glass exhibits non-Arrhenius behavior, and a sufficient model to represent such behavior requires at least three parameters. In this study, the major focus puts on three-parameter models only, able to describe the universal physics of the viscoelastic responses. A model in the form of Vogel-Fulcher-Tammann (VFT) equation is chosen given as follows: where T is the absolute temperature, and ℘ is each of the three temperature-dependent parameters of the Burgers model, ie, E 2 , η 1 and η 2 . In other words, each Burgers parameter can be determined at any temperature by the VFT function. Advantages of the proposed method are to increase the freedom of each model parameter when they are coupled with temperature, whereas all thermo-mechanical coupling parameters using the TRS property are restricted to a single temperature-dependent variable-the shift factor. Please note that the similar approach can be applied to any other rheological models instead of the Burgers model.
The VFT-fitting function in Equation (17) is chosen because each coefficient contains universal physics; though it is commonly known as an empirical estimation of viscosity with temperature dependence in the glass transition. In fact, the three coefficients, A, B, and T 0 , in the VFT form correspond, respectively, to (a) the viscosity at infinite temperature; (b) a function curve defining the fragility of a glass-forming material; and (c) a temperature of a glass-forming liquid at which its entropy, if the liquid is cooled sufficiently slow, coincides with the crystal entropy, defined as the Kauzmann temperature. 42 Understanding the physical meaning of each coefficient helps to minimize the number of fitting coefficients.

| Minimization of number of fitting coefficients
Using the VFT-fitting equation, the three temperaturedependent Burgers parameters are obtained, each of which contains three unknown coefficients. Thus, the total number of fitting coefficients is nine, including A(E 2 ), B(E 2 ), T 0 (E 2 ); A(η 1 ), B(η 1 ), T 0 (η 1 ) and A(η 2 ), B(η 2 ), T 0 (η 2 ). It is stressed that one can use all nine coefficients to determine each Burgers parameter as the temperature dependence.
For the benefit of computational effort, however, the next step is to find a less possible number of the fitting coefficients. First, all three Burgers parameters would have the same Kauzmann temperature since it is the unique property of a glass material; hence, T 0 (E 2 ) = T 0 (η 1 ) = T 0 (η 2 ) = T 0 . Besides, as observed in Figure 10, the variation of temperature-dependent viscosity η 1 (T) is mostly identical to that of E 2 (T). From the mathematical view using the VFT equation, both coefficients would have a similar coefficient B-the fragility, meaning that B(η 1 ) = B(E 2 ). Finally, we assume that the Kelvin component of the Burgers model has the same infinite viscosity, permitting A(η 2 ) = A(E 2 ). According to this thought, the number of fitting coefficient reduces to five, which are A(η 1 ), B(η 1 ), A(η 2 ), B(η 2 ), and T 0 .
As illustrated in Figure 10, the VFT model using a set of five fitting coefficients ℜ = ℜ A 1 , B 1 , A 2 , B 2 , T 0 yields excellent fits for all three Burgers parameters over the temperature range of investigation. It means that each Burgers parameter can be determined at any specific temperature. By using these parameters, eventually, either creep or relaxation responses can be directly identified without awareness of the relaxation functions and shift factors.

| Thermo-viscoelastic modeling procedure
Based on the discussions in Sections 4.1 and 4.2, a thermoviscoelastic modeling approach is established as in the following: The procedure starts with creep experiments conducted at a certain number of chosen temperatures n (eg, n = 4, Figure 3B) in the range of glass molding. Corresponding mold displacements u exp k over entire recording time t are acquired, while k is the index of the displacement (k = 1 − n). Goal of this study is the development of a computational method to estimate the creep displacements called u cal k , using the single set of Burgers parameters ℘. Objective of the parameter estimation procedure is to minimize the difference between the estimated and experiment displacements. More specifically, the estimation procedure is converted into an optimization problem so as to determine the unique parameter set ℘ by matching the estimated and experiment displacements simultaneously at all temperatures.
In fact, providing that the parameter set ℘ is known, the creep displacements are achieved according to the computational procedure presented in Section 3.2. Furthermore, the single set of Burgers parameters is able to compute all n displacements at different temperatures by using the temperature-dependent relation defined by Equation (17). Therefore, the computation of creep displacements over the experiment temperature range eventually requires determining the set ℜ of 5 unknown coefficients as discussed in Section 4.2.
A computational solver is developed for this purpose. Inputs fed into the solver consist of five unknown coefficients, the glass property (E, ) and the experiment data, ie, recording forces and displacements. Outputs are the computed displacements u cal k . At first, u cal k are computed using an arbitrary set of coefficients. Following, they are compared with the experiment data u exp k . Providing the difference is bigger than an acceptable tolerance, a new coefficient set ℜ is updated to correct the computed displacements iteratively until the following relation is satisfied: In this formulation, u exp k ( ) and u cal k ( ) are vectors containing the experiment and computed displacements, respectively, taken at a given time (0 < < t) and at the k-th temperature; t is the recording time of the displacement. The time-displacement data obtained from creep experiments are highly nonlinear as shown in Figure 2. The overdetermined system of linear equations cannot be solved for Equation (18). Therefore, the determination of the coefficient set ℜ is necessarily required a powerful nonlinear optimization. In this study, the Levenberg-Marquardt algorithm is utilized for the parameter estimation. Accordingly, the optimization procedure serves to estimate the computational displacements iteratively until the difference can be sufficiently negligible to satisfy Equation (18).
To calculate the displacements from the coefficient set, the solver consists of 3 modules. The first module is to determine the Burgers parameters from the coefficient set. At this point, note that the Burgers parameters can be converted into relaxation parameters in the form of Prony series, for the ease of model implementation into an FEM commercial software, such as ABAQUS. In the following step, as the Burgers parameters are known, the creep compliances are computed in the second module. Lastly, the third module is employed to calculate the displacements from the computed creep compliances. Detail of the computational procedure is described in the flowchart provided in Appendix B.

| Numerical validations
Numerical validations of the proposed method are conducted in ABAQUS. The FE model shown in Figure 5 is used to simulate the creep displacements at four testing temperatures. Viscoelastic parameters are provided in Table 3. The viscoelastic parameters are implemented for the simulation in ABAQUS via a UMAT user-subroutine. Figure 11A presents the simulation results, which greatly agree with the experiment data. Compared to the large discrepancies observed in Figure 8, the proposed method obviously enhances the simulation accuracy.
In order to demonstrate the efficiency of this method, further simulations are performed to predict the creep displacements at other temperatures that are not involved in the VFT-fitting function. For this reason, creep experiments from 650°C to 790°C with an increment of 10°C are validated. Using the same viscoelastic parameters in Table 3, FEM simulations are run at every experiment temperature. Figure 11B presents the simulation results and experimental data. Exceptionally good agreement is found within the entire temperature range chosen for the validation.
Furthermore, to intensify the appropriateness of the method, other glass types, commonly used in PGM, are investigated. Beside the borosilicate glass SUPRAX ® 8488, two other glass categories of SCHOTT, aluminosilicate Xensation ® 3D and infrared chalcogenide IRG 27 (As 2 S 3 ), are chosen. Thermomechanical properties of these glass materials are provided in Table 1. Figures 12 and 13 introduce the creep experiments and numerical validations of the Xensation ® 3D and IRG 27, respectively. For those plots, the dashed lines in Figures 12A and  13A present the numerical verifications at the testing temperatures chosen for obtaining the Burgers parameters by the proposed method. The VFT-fitting functions of those parameters are shown in Figures 12B and 13B, and the corresponding values are provided in Table 3. Besides, the solid lines exhibit the validation of the creep displacements at temperatures that are not involved in the fitting functions. These results demonstrate that accuracy of the numerical predictions at other temperatures within the glass molding range is satisfied. The validation results of three different glass categories strengthen the validity of the proposed method in characterizing and modeling the thermo-viscoelastic nature of glass for FEM simulation in the glass molding process. In addition to considering the VFT model, it is emphasized that other existing equilibrium viscosity models can be valid for the proposed method. To demonstrate, the MYEGA viscosity model, recently introduced by Mauro et al 43 is investigated. Similar as the VFT form, this model contains three fitting parameters that are used to obtain the temperature-dependent Burgers parameters. Figure 14 compares the simulation of creep displacements of the MYEGA and the VFT models where exceptionally identical results are found. In future work, the validity of the equilibrium viscosity models needs to further examine a broader temperature range, particularly for temperatures near T g and below.
Finally, the methodology is independent from the choice of material models defined by simulation users. For example, the generalized Maxwell model with multiple Maxwell terms n ≥ 3 are commonly used in many previous studies. We emphasize that the same modeling procedure can be truly adopted for other user-defined material models.

RESEARCH
An accurate model to characterize the thermo-viscoelastic properties is among the most successful keys for the simulation of the PGM process. In this paper, present restrictions on the modeling of the viscoelastic material behavior have been resolved. For the simulation of the viscoelastic responses, a correct implementation of the relaxation parameters derived from the creep experiment data is elucidated. The finding gets rid of the mismatch between the simulation and experiment data of the creep displacements. For modeling the temperature-dependent viscoelasticity, a detailed examination reveals the deficiency of the TRS assumption over a wide temperature range of glass molding. To overcome this deficit, the paper presents a novel method for modeling the thermoviscoelastic material behavior of glass. The proposed method permits a determination of the thermo-viscoelastic parameters without using the relaxation functions and shift factors. The constitutive Burgers model demonstrates its sufficient representations of the viscoelastic responses over the range of molding temperature chosen in this study. Temperature is incorporated into each Burgers parameters in the form of VFT function to account for the temperature-dependent viscoelasticity. The method illustrates that a single set of five coefficients satisfactorily characterizes the thermo-viscoelastic properties of the glass. By multiple curve fitting of the creep displacement data, the parameter set is achieved. Using this single set, the viscoelastic responses, either creep or stress

MATHEMATICAL EQUIVALENCE BET WEEN THE BURGERS MODEL AND THE T WO-TERM PRONY SERIES
Constitutive equation of the Burgers model in the differential form is given as (Findley et al 29 ): Using Laplace transformation, the solution to creep and stress relaxation responses can be obtained: and where q 1 ,q 2, and A are the parameters defined as follows: with and 1 , 2 are the relaxation times, given as follows: Following, shear parameters of the Burgers model are calculated by using Lamé constants: