Rheological measurement of the nonlinear viscoelasticity of the ABS polymer and numerical simulation of thermoforming process

The thickness distribution of thermoformed products is greatly affected by the viscoelastic behavior of the extruded polymer sheet. In this work, linear and nonlinear rheological experiments are carried out to characterize the viscoelastic properties of acrylonitrile-butadiene-styrene sheets under thermoforming conditions including a wide range of temperatures, strains, and strain rates. First, aspects of linear viscoelasticity such as the storage modulus and loss modulus are measured by small-amplitude oscillatory shear experiments. The discrete relaxation spectra and the Williams-Landel-Ferry parameters are obtained from the constructed linear master curves. Then, nonlinear time-dependent extensional viscosity is measured by uniaxial extensional experiments. The parameters of the damping function are evaluated using an optimization method. In addition, the effect of the orientation of the polymer is analyzed. The uniaxial extensional stress and viscosity in the extruder direction demonstrate higher resistance against tearing and extreme thickness reduction during processing. Finally, the linear and nonlinear input parameters for the numerical simulation are prepared. Numerical simulations are performed using the Wagner model with the obtained nonlinear viscoelasticity. The thickness distribution in thermoformed ABS sheets, obtained numerically, shows good agreement with the experimentally obtained values.


Thermoforming
Thermoforming is a manufacturing process used to convert polymer sheets into thin-walled plastic products with desired shapes. In the thermoforming process, the polymer sheet is deformed by the heat and pressure applied. The sheet, heated above the glass transition temperature, has both viscous and elastic components. The sheet undergoes rapid and large deformation during the process, and this results in a change in viscoelasticity, which greatly affects the thickness distribution of the thermoformed product. Therefore, it is very important to investigate the viscoelastic behavior of polymer sheets under thermoforming process conditions such as strain, strain rate, and temperature.
The advantages of thermoforming are low cost in mold design, precise shape, low pressure, automation for mass production, and the ability to use a variety of applicable polymer materials. Thermoforming applications are classified into thin-gauge and heavy gauge types, depending on the thickness of the product. The former group includes packaging containers with thicknesses of less than 1 mm, such as food, beverage, and pharmaceutical and semiconductor containers. The latter group includes products with thicknesses from 1 to 10 mm and includes automotive, ship, and aircraft parts that require high durability. Polymer material undergoes large and rapid deformation during the thermoforming process, and this can result in nonuniform thickness distribution, excessive thinning and a decrease in mechanical properties. Typically, nonuniform thickness distribution occurs more frequently in heavy gauge thermoforming, which involves deep drawing molds with complex shapes, as compared with thin gauge thermoforming.
The thermoforming process consists of three main steps: heating, pre-stretching, and forming, as shown in Fig. 1. In the heating step, shown in Fig. 1a, the polymer sheet is heated above the glass transition temperature (T g ) by infrared radiation, which causes it to lose its mechanical strength; eventually, it sags due to gravity. In the pre-stretching step shown in Fig.  1b, the polymer sheet is pre-stretched to improve the uniformity of thickness distribution in the thermoformed product. The edge area of the polymer sheet is fixed by a clamping tool and vacuum pressure is applied to inflate the sheet surface like a bubble. In the forming step shown in Fig. 1c, the polymer sheet undergoes extensional deformation due to contact with the mold, and local cooling occurs where areas of the polymer sheet are in contact with the lower temperature mold. This reduces the mobility of the styrene-acrylonitrile (SAN) matrix, and the extension of the polymer sheet stops. The deformation behavior that occurs in each step of thermoforming, such as sagging, expansion, extension, and cooling by the mold, affects the thickness distribution of the thermoformed product.
In the thermoforming process, polymer materials exhibit deformation-and temperature-dependent nonlinear viscoelasticity depending on time and space. Under particular process conditions, the internal stress grows or relaxes over time. The characteristics of the polymer material are affected by temperature. In the glassy state, the polymer material is brittle and only small deformation is possible. However, since the material is ductile in the rubbery state above the T g , large deformation is possible. In addition, stress can be increased or decreased depending on the strain rates applied to the material. Polymer materials can have different mechanical properties depending on their molecular weight and structure, and polymer sheets produced by the extrusion process can exhibit strong anisotropy due to resulting molecular orientation. The thickness distribution of thermoformed products is greatly affected by the viscoelastic characteristics of the polymer materials.
Industrial thermoforming processes are typically designed to improve quality and uniform thickness distribution, increase production efficiency, and reduce material costs. These objectives can be affected by various parameters of the thermoforming process, such as vacuum pressure, sheet temperature, mold temperature, ambient temperature, convective heat transfer coefficient, mold speed, mold material, coefficient of friction, and thermal conductivity. Trial-and-error-based experiments have been conducted to investigate the effects of process parameters. However, that approach requires significant effort and cost, and it can be difficult to analyze the deformation behavior of the polymer materials in detail because there are so many parameters. In addition, a lack of detailed theoretical background can lead to problems, such as unoptimized process, process instability, material waste, and low-quality products. A numerical simulation is an excellent tool for conducting in-depth analyses of the deformation of polymer materials and the thickness distribution of thermoformed products. To perform numerical simulations, a nonlinear viscoelastic constitutive model is required, with material properties under process conditions.

Literature review
Many previous works have been devoted to the rheological measurement of linear and nonlinear viscoelastic properties, as well as numerical simulations with constitutive models. Wang  [1] performed uniaxial tensile experiments using a Meissner-type rheometer to measure the properties of ABS and proposed a material model that considers the effects of strain rate, strain hardening, and temperature sensitivity. Lau et al. [2] showed that ABS exhibited high resistance to sheet sag because of its high viscosity and high modulus of elasticity. Because of its broader molecular weight distribution, polypropylene had a sag resistance lower than that of ABS. Aoki et al. [3] studied the effect of the dispersion of rubber particles in ABS polymer on damping function, which indicated stress reduction with material deformation. It was found that the damping function became stronger with increased rubber particle content in the ABS sample containing well-dispersed rubber particles. Lee et al. [4] measured the rheological properties of ABS, such as elongational viscosity, strain hardening, and strain softening using a Meissner-type extensional rheometer and mechanical spectrometer. They reported that the strain-hardening behavior of the polymer material was important for uniform thickness distribution. Münstedt et al. [5] showed that the strain hardening of polymer materials, such as polyethylene and polypropylene, had a significant impact on the uniformity of deformation. The innovative thermoforming process has been recently developed to form composite sheets by means of polymer melt pressure during the injection phase [6]. This process is an integration of the thermoforming and injection molding processes and can result in significant reductions of tooling, machinery, and operational cost [7]. In general, the extruded polymer sheets used in thermoforming exhibit anisotropic stresses along the machine direction or in the transverse direction. Chan and Lee [8] reported that the orientation stress of an extruded polypropylene sheet depended strongly on the draw ratio and that higher melt temperature reduces the orientation effect. Stephenson et al. [9] conducted an experimental study of the sag of an extruded polypropylene sheet and found that the extension ratio in the transverse direction was higher than in the machine direction. The polystyrene was found to be temperature-sensitive in the thermoforming range, and temperature effects should be included to obtain a better understanding of sheet sag.
Two different types of constitutive models can be considered to perform numerical simulations of the thermoforming process. Hyperelastic models are commonly used because of their simplicity, and viscoelastic models are employed to describe the deformation-and temperature-dependent nonlinear viscoelasticity of polymer materials. Koziey et al. [10] compared the hyperelastic Ogden model and the Kaye-Bernstein-Kearsley-Zapas (K-BKZ) viscoelastic model in a numerical simulation of thermoforming. The Ogden model provided reliable results for molds with shallow depths and simple shapes, while the K-BKZ model showed good agreement with thickness distribution using deep drawing molds with complex shapes. Novotný et al. [11] reported that the hyperelastic model was not suitable to describe the deformation of polymer materials in thermoforming in detail and that the nonlinear viscoelastic model was capable of accurately representing the behavior of polymer materials. The K-BKZ model was used to describe viscoelasticity for large deformations. Sala et al. [12] performed numerical simulations with the G'sell viscoelastic model of high impact polystyrene (HIPS) thermoforming for a refrigerator liner. A 20% discrepancy between experimental measurements and the numerical results was found, which it may be possible to improve by integrating the heat transfer and the material behaviors. Dong et al. [13] used hyperelastic models such as the Mooney-Rivlin and Ogden models to carry out numerical simulations of bubble inflation of a polymethylmethacrylate sheet. Anisotropic material properties can have a significant impact on simulation results when applying an isotropic material model, but there were no significant differences in PMMA ratings after shrinkage testing.
There have been several studies that have demonstrated the influences of various process parameters, such as heating, temperature of sheet/mold/plug, speed and shape of mold/plug, pressure in air/vacuum, and coefficient of friction on the thickness distribution of thermoformed products. Ayhan and Zhang [14] studied the effects of sheet temperature, air pressure, and heating time for simple-shaped food containers in plug-assisted thermoforming using multilayered polymer sheets. The most important parameter affecting the wall thickness distribution was sheet temperature. Bhattacharyya et al. [15] used the grid strain analysis technique to quantify the difference in the strain distribution during the thermoforming of composite sheets made of wood fiber and polypropylene. They reported that the forming temperature and the size of the clamp frame had the most significant effect on thermoformability. Chen et al. [16] investigated the thickness distribution of thermoformed thin films using 0.125 mm and 0.2 mm polycarbonate. They found that thickness at the sidewalls increased with increasing mold temperature, preheating temperature, plug depth, and holding time but decreased with increasing plug speed. McCool and Martin [17] studied plug-assisted thermoforming for a simple shaped cup using a 1.25-mm thick HIPS polymer sheet. Plug displacement, sheet temperature, plug temperature, and plug shape greatly influenced thickness distribution, but plug speed and air pressure had relatively small influence.
Collins et al. [18] studied the role of contact between the sheet and several plug materials such as syntactic foam, polyoxymethylene, and aluminum. The experimental results showed that under T g the coefficient of friction (COF) is in the range of 0.1 to 0.3, but COF increased sharply with increasing temperature. At a thermoforming temperature above T g , the sheet material became sticky and could not slip on the plug surface. O'Connor et al. [19] found that in the plug-assisted thermoforming of polypropylene cups the contact friction between the sheet and plug were the most sensitive parameters to control thickness distribution. Morales et al. [20] measured the coefficient of friction between a HIPS sheet and plug materials such as velocities of aluminum and steel in thermoforming. The COF increased with increasing sheet temperature, but there was no significant change with plug speed. Marathe et al. [21] studied the temperature dependence of COF in plug-assisted thermoforming using 3.0 mm thick HIPS. It was found that the COF rapidly increased when the temperature was higher than the T g , resulting in a significant reduction in slip and increasing uniaxial stretching at the clamp.
Furthermore, many studies have highlighted the significance of the heating step in the thermoforming process, because the thickness distribution is affected by the temperature-dependent viscoelastic properties of the polymer materials. Yousefi et al. [22] focused on accurate measurement of heat transfer parameters such as heat transfer coefficient, atmospheric temperature, emissivity, view factor, specific heat, and thermal conductivity for the 1.6-mm thick ABS sheet. Duarte and Covas [23] attempted to solve the inverse heating problem in thermoforming, minimize the temperature gradients across the sheet thickness, and obtain uniform thickness distribution. Because of significant sheet deformation due to sheet sag during the heating step, Wiesche [24] used a dynamic view factor to perform radiative and conductive heat transfer simulations. Labeas et al. [25] studied the optimization of the heating step by conducting a simulation involving heat transfer between thermoplastic carbon/polyetherimide composite sheets. Various parameters, such as the number of heating elements, and location, and heat flux of infrared heating lamps were investigated for optimal heating. Giacomin et al. [26] reported that the complicated processability of thermoforming comes from temperature gradients, either across the sheet or through its thickness, due to sheet sag. Ghobadnam et al. [27] found that different heating times are required depending on the shape of the mold. Increasing the heating time was important to achieve a more uniform thickness distribution for the female mold. However, for male molds, the heating time had to be reduced.
From the previous literature, it is clear that the thickness distribution of thermoformed products is dependent on the behavior of the polymer materials. In addition, the investigation of polymer materials should be conducted under process conditions including a wide range of temperatures, strains, and strain rates. Moreover, numerical simulation is essential to a better understanding of the thermoforming process via the study of the effects of various process parameters. Hence, rheological measurement of the nonlinear viscoelasticity of the polymer material should be considered.

Objectives
The objective of this study is to investigate the nonlinear viscoelastic characteristics of the extruded polymer sheet used in heavy gauge thermoforming. Rheological measurements were carried out under thermoforming conditions including a wide range of strains, strain rates, and temperatures. The linear and nonlinear viscoelastic properties of the extruded ABS sheet were extensively investigated, including by the constitutive model, and the effects of orientation in terms of viscosity, stress, and homogeneous deformation were analyzed. Numerical simulations based on the finite element method were performed with the obtained nonlinear viscoelastic parameters. The numerical results were validated by comparing them with the thickness distribution of products thermoformed by an industrial thermoforming machine. This work provides important insight into the deformation behavior of polymer materials that undergo large and rapid deformations in thermoforming.

Nonlinear viscoelastic constitutive model
Generally, polymers are viscoelastic materials that dissipate a portion of the applied energy and simultaneously remember the deformation history. Polymer materials exhibit linear viscoelasticity in small deformations, but gradually display nonlinear viscoelastic behavior under large strain and high strain rates over a wide range of temperatures. Considering that deformations are large and rapid in most processing operations, it is necessary to describe the nonlinear behavior of polymers accurately [28]. To this end, many constitutive equations have been proposed in integral or differential forms [29]. Among the constitutive equations, the K-BKZ equation, developed by Kaye [30] and Bernstein et al. [31], has been extensively evaluated due to its wide applicability to polymer materials. In particular, it was reported that the K-BKZ model results in a credible agreement with the experimental uniaxial extensional behavior of ABS [4]. Thus, in this work, the K-BKZ model was selected to investigate the linear and nonlinear behavior of ABS.
The K-BKZ equation is a class of equations because its predictability depends on the choice of the kernel function. A simple and useful equation form was suggested by Wagner, which is called the Wagner model [32]. The Wagner model is expressed as follows: where σ is an extra molecular stress tensor, M is a kernel function that depends on both time and strain, B is a Finger strain tensor, and I 1 and I 2 are the first and second invariants of the Finger strain tensor, respectively. Wagner proposed the kernel function as the product of a time-dependent memory function and a strain-dependent damping function, which is called time-strain separability.
where m(t − t′) is the time-dependent memory function used to explain the linear viscoelasticity and h(I 1 , I 2 ) is the straindependent damping function used to describe the nonlinear viscoelasticity. The time-dependent memory function and strain-dependent damping function should be specified in this model to explain the linear and nonlinear viscoelastic behaviors.
The memory function is presented as the sum of the discrete relaxation modes with discrete relaxation moduli g i and discrete relaxation time λ i , Generally, N sets of relaxation moduli are called the discrete relaxation spectra. The values g i and λ i can be experimentally obtained by fitting the following equation to small amplitude oscillatory shear (SAOS) data, Various damping functions have been suggested based on theoretical and empirical approaches [33][34][35][36][37][38]. The equation forms for damping functions are summarized well in a recent review by Rolón-Garrido and Wagner [39]. The damping function h(I 1 , I 2 ) as a function of the invariants of the Finger strain tensor has a value between 0 and 1.
When h(I 1 , I 2 ) is one, the materials exhibit linear viscoelasticity. In this case, Eq. (1) reduces to the Lodge model developed by Lodge [40]. As h(I 1 , I 2 ) goes to zero, the nonlinear viscoelasticity of the materials becomes noticeable. The damping function of the Wagner model, used in this study, is expressed as where α and β are the nonlinear parameters of the damping function, obtained from the uniaxial extensional experiments. The first and second invariants of the Finger strain tensor for uniaxial extensional flow are as follows: where ε h is the Hencky strain, defined as the product of extensional rate and time (≡ ε t). By combining Eq. (1) with Eqs.
3 Rheological measurement of material properties

ABS polymer
ABS is a copolymer composed of acrylonitrile, butadiene, and styrene. It is widely used in thermoforming because it has high impact resistance, rigidity, gloss, hardness, and temperature stability. Furthermore, ABS has high melt elasticity because butadiene rubber particles of several, hundred nanometer sizes are dispersed in the SAN matrix. The commercial ABS copolymer used was provided by LG Chemistry Corporation (South Korea). The SAN matrix (M w = 206,500 g/mol) consisted of 31.8 wt% acrylonitrile and 68.2 wt% styrene. SAN-grafted butadiene rubber phase of 10 wt% was added to the matrix. The rubber particles have a narrow size distribution with an average diameter of 230 nm.

Glass transition temperature
Glass transition temperature (T g ) can be used to determine the thermoformable temperature range. ABS is an amorphous polymer and, thus, does not have a clear melting point. ABS is hard and brittle in the glassy state below T g . In the rubbery state above T g , the mobility of the polymer chain increases, and thus ductility increases and rigidity decreases. Accordingly, in this work, the thermoforming temperature range was set from T g + 30 to T g + 60°C.
To measure the T g of the ABS used, a dynamic mechanical analysis (DMA) experiment in cantilever mode was performed using a DMA Q800 (TA Instruments). Specimens were obtained by cutting the extruded ABS sheet into dimensions of 35 mm long, 13 mm wide, and 3.4 mm thick. The elastic and viscous responses (E′ and E″) were measured at three frequencies of 1, 5, and 10 Hz in a range of temperatures from 30 to 150°C, with a heating rate of 2°C/min. From Fig. 2, it is clear that the viscoelasticity of the ABS has temperature dependence. The E′ of the ABS did not change significantly up to 100°C, while E″ changed sensitively with increasing temperature. T g was obtained at the peak of the tan δ Ε (≡ E″ / E′) curve. Near T g , the polymer material exhibited a transition from glassy to rubbery state, and thus E′ and E″ changed dramatically. The measured values of T g at the three frequencies are summarized in Table 1.

Linear rheological measurements
At low strain and strain rates, a polymer material exhibits linear viscoelastic behavior. However, as strain and strain rates increase, the material response deviates from linear viscoelasticity and gradually moves to nonlinear viscoelastic behavior. SAOS experiments were performed to characterize the linear viscoelastic properties of ABS over a wide range of temperatures, from 130 to 220°C, with an angular frequency ω of 10 −1 to 10 2 rad/s, and at small shear strain amplitude, γ 0 = 0.5%, corresponding to the linear viscoelastic condition.
The linear viscoelastic storage modulus (G′) and loss modulus (G″) are measured in SAOS tests. When an input strain signal is a simple sine function, the shear stress signal under linear oscillatory shear is expressed as follows: where G′(ω) is the storage modulus, G″(ω) is the loss modulus, and tan δ is the ratio of G″ (ω) to G′(ω). SAOS experiments were carried out using a rotational rheometer (ARES-G2, TA instruments) with parallel-plate geometry and a diameter of 13 mm. A circular specimen with a thickness of 3.4 mm and a diameter of 13 mm was cut from the extruded ABS sheet. Each experimental test was repeated three to five times, and one representative test was selected for the analysis. Based on the time-temperature superposition principle, the experimental data at different temperatures were shifted to construct a linear master curve. The data were shifted along the frequency axis, and the horizontal shift factors a T were fitted with the Williams-Landel-Ferry (WLF) equation.
where T is the current temperature, T r is the reference temperature, and C 1 and C 2 are the parameters of the WLF model. The parameters C 1 = 4.67 and C 2 = 105.14 K were obtained at a T r of 170°C. As depicted in Fig. 3, the shift factors are described well by the WLF equation. The obtained linear master curve is shown in Fig. 4, with the prediction from Eq. (4). Ten sets of the discrete relaxation time λ i and discrete relaxation modulus g i were obtained and are summarized in Table 2.

Nonlinear rheological measurements
Linear viscoelastic behavior can be analyzed using G′(ω) and G″(ω) when the deformation is small. However, the linear viscoelastic theory is no longer valid when polymer materials undergo large and rapid deformation. As strain increases, material response gradually departs from the linear behavior, and nonlinear behavior such as strain hardening and strain softening becomes noticeable. It is therefore very important in the thermoforming process to analyze the nonlinear viscoelastic behavior of polymer materials.
To this end, the uniaxial extensional viscosity of ABS was measured using an ARES-G2 equipped with an extensional viscosity fixture. Figure 5 shows the principle of uniaxial extensional viscosity measurement. One end of the specimen is connected to a fixed drum where the force is measured. The other end is connected to a rotating drum pulling the specimen. The rotating drum moves along a circular orbit around the stationary drum to apply uniaxial stretch with a constant extensional rate to the specimen.
Since the polymer sheet used in thermoforming is extruded, a strong orientation occurs inside the sheet with respect to the extrusion direction, which may exhibit significant anisotropy. The oriented polymer sheet has a high resistance to deformation in the extrusion direction and exhibits a weak resistance in the transverse direction, which may have a large influence on thermoforming. To investigate the orientation where σ E t; ε ð Þ is the extensional stress in terms of time and strain rates, F(t) is the force, and l 0 and w 0 are the length and width of the specimens, respectively. To check the accuracy of the extensional measurements, η þ E t; ε ð Þ was compared with the linear viscoelastic behavior, calculated using λ i and g i from the SAOS experiments, Figure 6 shows the uniaxial extensional viscosity curves of the extruded ABS sheet for different extensional rates at various temperatures. The dashed line indicates the linear viscoelastic response of η þ E t; ε ð Þ (Eq. 12). The filled symbol (type A) and open symbols (type B) indicate the experimentally obtained η þ E t; ε ð Þ. In the initial stage, η þ E t; ε ð Þ of both types showed good agreement with the linear viscoelastic response when the Hencky strain was small. As the strain increased, η þ E t; ε ð Þ gradually deviated from η þ E t ð Þ, with significant strain hardening. Strain-hardening behavior became more remarkable as the extensional rates increased at a fixed temperature. Accordingly, rapid deformation can have a positive effect on the thickness distribution of thermoformed products, by preventing excessive thinning. As the strain increased further, softening behavior was observed, indicating that melt failure occurred near the maximum strain of 3.4. These behaviors indicate that a nonlinear constitutive model is required to describe the nonlinear viscoelasticity of the ABS under thermoforming conditions, such as large strain and/or high strain rate.
As the temperature increased, it was found through analysis that the increase in mobility of the SAN matrix influences the uniaxial extensional behavior of ABS sheets. The η þ E t; ε ð Þ curves at 160°C displayed strong strain hardening behavior. However, as the temperature decreased, the strength of strain hardening also decreased. The η þ E t; ε ð Þat 190°C even showed strain-softening behavior and almost linear viscous characteristics at 0.1 s −1 . At high temperatures, low η þ E t; ε ð Þ and strain softening can lead to unstable deformation, such as a bubble explosion in the pre-stretching step and tearing in the forming step.
The stress-strain curves were plotted to clearly show the effect of orientation of the ABS sheet, as shown in Fig. 7. Within for students the experimental conditions, type A mostly exhibited a higher stress value and earlier strain hardening than type B. This difference between the two types has a significant effect on the thickness distribution of the thermoformed products. The earlier strain hardening of type A provides high resistance to the applied deformation, preventing extreme thickness reduction or tearing due to high frictional force. In contrast, type B has a lower resistance to deformation at ε h < 3. Type B is also expected to display enough strain hardening at ε h > 4. However, considering that the maximum ε h value during the thermoforming process is less than 4, as determined by Throne [41], a more uniform thickness distribution can be achieved when the ABS chains are oriented toward the extension direction. Extensional viscosity, strain hardening, and the effects of orientation on stress  and strain can be very useful information in polymer processing.
The strain hardening behavior of ABS can be controlled by the amount and hardness of butadiene particles. Takahashi et al. [42] investigated the effect of the hardness of butadiene particles for a series of ABS samples that had similar particle size and dispersion state. Notably, the strain hardening became weaker with increasing content of soft or hard particles. These results suggest that the degree of strain hardening is strongly affected by the particle hardness. Thus, the role of butadiene particles on strain hardening behavior is highly important. Compared with their results, a particle content of 10 wt% is  Fig. 5 Schematic of the extensional viscosity fixture to measure the uniaxial extensional viscosity adequate for thermoforming processes because it generates enough strain hardening to prevent excess thickness thinning.

Orientation effects of extruded sheet
Strain hardening of polymer material has a positive effect on homogeneous deformation in the polymer process, and this was quantitatively analyzed using the Considère criterion [43].
The Considère criterion has been used to predict the failure of polymer melts and solutions. Lee et al. [4] in particular used the Considère criterion as an effective indicator of the nonuniform deformation of ABS in uniaxial extension. Figure 8 shows the critical strain (ε c ) obtained with Eq. (13) along with force and strain curves for type A and type B, respectively. The ε c values are summarized in Table 3. In all the results, ε c increased with increasing extensional rates at a lower temperature.
The ε c value represents the maximum Hencky strain at which samples are deformed homogeneously. Thus, it can be regarded as a measure to determine the operating limit in polymer processes such as thermoforming, which are important for achieving uniform thickness distribution. At 180 and 190°C, the two types displayed similar ε c values (Table 3). However, with decreasing temperature, the ε c value of type A shifted to two, while type B remained less than one. Thus, type A can undergo homogeneous deformation until higher ε c , compared with type B. From these results, it can be seen that a critical strain of an extruded ABS sheet exists at which homogeneous deformation can be maintained using strain, strain rate, and temperature.
Taken together, the thickness and uniformity of the thermoformed products highly depend on the extensional rheology of the polymer material. The different extension directions with respect to chain orientation, in particular, can result in different nonlinear responses, which leads to the final Fig. 7 a-d Stress-strain curves of ABS for various extensional rates in the temperature range of 160 to 190°C. Under the same conditions, the stress of the type A specimen is larger than that of type B specimen properties of the thermoformed products. As a result, it was speculated that uniform thickness distribution and high tear resistance can be obtained when the sheet is pulled in the extrusion direction.

Parameters of damping function
The Wagner model from Section 2 was exploited to describe the uniaxial extensional behavior of the ABS. The specific form of the equation is as follows: where s = t − t′. The linear viscoelastic relaxation modulus G(t) and memory function m(s) were calculated from the discrete relaxation spectra (Eq. 3). Equation (6) was used for the damping function h(t) and h(s). To determine the nonlinear parameters α and β in the model, the following objective function (RSS) was optimized numerically.
where F exp i is the experimental data, F model i is the Wagner model predictions, and D is the total number of data points.
During the thermoforming process, because deformations occurred mainly in the machine direction of the extruded ABS sheet, the parameters α and β were obtained using the type A specimen. The reference temperature of 170°C was selected to describe the temperature-dependence, using the WLF equation. Figure 9 shows the uniaxial extensional viscosity curves of experiments and the predictions from the Wagner model using the obtained parameters. The Wagner model prediction for η þ E t; ε ð Þ can be used to quantitatively understand the characteristics of the experimental results. At a short time, the model prediction coincided with the linear viscoelastic response, validating the accuracy of the numerical solutions. As time increased, the Wagner model gave good predictions for strain hardening behavior, especially the onset time of strain hardening. In contrast to the Lodge model, which predicts infinite growth of the extensional viscosity (Lodge, 1956), the Wagner model allows one to obtain the maximum viscosity value as a function of the damping function term. However, at an extensional rate of 0.1 s −1 , the model predicted higher η þ E t; ε ð Þ, and stronger strain hardening than the experimental data. Accordingly, a complete description of the experimental data was not clear when using the parameter sets of α and β. The validity of the parameters will be discussed in Section 4.3; the results of numerical simulations using the Wagner model can reasonably explain the deformation of the extruded ABS sheet compared with that in the experimental results.

Numerical setup
Three-dimensional simulations based on the finite element method were performed using T-SIM (version 4.83, Accuform), a commercial code package for thermoforming. The deformation of the polymer materials and the thickness distribution of the thermoformed products were numerically investigated. Numerical simulations require a nonlinear viscoelastic constitutive model that can describe the relationship between stress and deformation. The discrete relaxation spectra, the WLF equation, and the Wagner model were utilized in the numerical simulations.
A mold for a refrigerator inner case using heavy gauge and deep drawing thermoforming was designed using the commercial CAD software NX (version 8.0, Siemens). This refrigerator case had a complex shape because it contained structures to increase the structural rigidity of the final product and to mount various parts. Figure 10 shows the shape and dimensions of the mold. The width, depth, and height of the mold were 840 mm, 740 mm, and 1620 mm, respectively. The mold consisted of an upper compartment 820 mm in length and a lower compartment of 760 mm length. The partition of the narrow gap between the compartments prevents the selfcontact of the polymer sheet during bubble inflation in the prestretching step. The frame clamps the polymer sheet during the entire process. Figure 11 shows the finite element model used for the numerical simulation. The polymer sheet was approximated as triangular membrane elements because the thickness of the sheet was very thin compared with the length or width, and thus the velocity and stress changes in the thickness direction could be ignored. Approximately 100,000 finite elements were established by performing a grid dependency test; test results showed excellent numerical convergence. Because the mold was symmetrical, only half of the polymer sheet was used for efficient calculations. Tools such as molds, partitions, and frames are considered rigid bodies that are not deformed and penetrated. The boundary conditions for the sheet were set to have zero velocity for the clamped edge and a symmetrical condition for the half-line. The initial temperature fields for  Figure 12 shows the process conditions such as differential pressure and mold position over time. From 0 to 4 s in the prestretching step, differential pressure was applied to inflate the polymer sheet. The mold stayed at the initial position until t = 2 s and then reached the maximum position at t = 4 s. From t = 4 s to t = 8 s in the forming step, the polymer sheet was stretched by contact with the mold; then, a vacuum was applied to the sheet through fine holes in the mold to obtain the desired shaped products.
During the process, the mold temperature can be increased by repeated contact with the heated polymer sheet, and so water flow channels were fabricated in the mold for cooling. The heat transfer between the sheet and the mold is expressed as follows.
where q is the heat flux, U is the overall heat transfer coefficient, T sheet is the temperature of the sheet, and T mold is the temperature of the mold.

Numerical results
The numerical results at each step of the thermoforming process are illustrated in Fig. 13. The polymer sheet was assumed to be preheated and sagged due to gravity, as shown in Fig. 13a. Next, vacuum pressure was applied to inflate the polymer sheet like a bubble for more uniform thickness distribution, as shown in Fig. 13b. This helps to distribute the material in contact with the rear surface of the mold into other areas before the forming step. As plotted in Fig. 13c, the inflated sheet begins to touch the edge of the mold. The small contact area between the sheet and the mold causes partial slips, not complete sticking. Then, the sheet continues to experience extensional deformation due to the movement of the mold. Finally, the mold stops at its predesigned position, and then vacuum pressure is applied through fine holes in the mold to fully contact the polymer sheet with the entire surface of the mold. A thermoformed product with the desired shape is obtained; the thinnest area was found near the edge of the mold. Figure 14 shows the thickness distribution on the symmetric line A-A' of products obtained from the experiments and numerical simulations. Thermoforming products were then manufactured using industrial thermoforming machines with a deep drawing mold with a complex shape. The same experimental and numerical conditions were used for comparison, as described in Section 4. were free surfaces in contact with the mold. They were mainly uniaxially deformed, and significant thickness reduction was observed because these regions continued to deform until the mold stopped. In contrast, the rear region was mainly biaxially deformed and the thickness change was relatively small. The hot sheet was cooled by contact with the edge of the mold so that it became difficult to pull the material from the rear region to the top and bottom regions. Due to a combination of cooling and partial slips, nonuniform thickness distributions were found in both the experimental and numerical results. Maximum deviation between experimental data (1.1 mm) and numerical results (0.77 mm) was observed at the edge between the top region and rear region in the upper compartments because temperature-dependent COF was not included. Figure 15 provides a visualization of the deformation of the side walls: the experimental and numerical thickness distributions are compared. Before deformation, a grid was drawn on the sheet with a height of 62.5 mm and a width of 125 mm. After deformation, the height of the grid in the upper compartment was 208 mm in the experiments (Fig. 15a) and 170 mm in the numerical results (Fig. 15b), and the height in the lower compartment was 214 mm in the experiments (Fig. 15a) and 158 mm in the numerical results (Fig. 15b). High extensional deformation at the side walls was observed in both results, and the experimental results showed a larger change than that of the numerical results. As a result, in the sidewalls, the experimental thickness was lower than the numerical results (500-1000 mm in line B-B' and line C-C', corresponding to the upper and lower compartments, respectively). This is due to the limitations of the constitutive model, which cannot describe anisotropy, such as the effect of the orientation.
As shown in Figs. 6, 7, and 8, the extruded polymer sheet was sensitive to the applied deformation along a specific direction. However, isotropic models, such as the Wagner model, used in this work, do not take this effect into account. Numerical simulations of thermoforming could be improved  by the development of anisotropic viscoelastic constitutive models that are capable of a prediction considering the orientation of the polymer chains. The thickness distributions in the numerical simulations showed good qualitative agreement with the experimental results, but some discrepancies remain. This can be caused by the use of a constant value of COF in the simulation. As shown in previous works [18][19][20][21], COF is a key parameter to control the thickness distribution in thermoforming and is a function of temperature. Since the inflation of the polymer sheet leads to a free surface in the pre-stretching step, the influence of the coefficient of friction on the deformation behavior is not large. On the other hand, the extensional deformation of the polymer sheet can be greatly affected by contact friction in the forming step. As Marathe et al. reported [21], COF is highly dependent on temperature and increases very steeply above T g . The measured COF of HIPS was low for T < 100°C, whereas it increased considerably at and above 100°C in their work. It is expected that the behavior of sheet sticking or sliding on the surface of the thermoforming tools can significantly change with temperature change. For example, using a high COF can increase the thickness, as the sheet sticks to the surface of the mold, and deformation stops in that area. Contact friction is the most important factor affecting the thickness distribution of thermoformed products. Hence, for a more accurate prediction of the thickness distribution, temperature-dependent COF of the polymer sheet should be incorporated into the numerical simulation.
Furthermore, the assumption of a uniform temperature field at the initial state for each compartment can cause the differences between the experimental results and the simulated predictions. In the heating step, the industrial thermoforming machine uses hundreds of heating elements to heat the polymer sheet, and it is difficult to obtain a uniform temperature field of the sheet. Small temperature differences across the sheet surface can exist, and temperature-dependent viscoelastic properties of the sheet can lead to error.

Conclusion
In this work, using a deep drawing mold with a complex shape, rheological measurements were carried out to characterize the linear and nonlinear viscoelasticity of the ABS used in heavy gauge thermoforming. Moreover, the effect of the orientation of the extruded polymeric sheets was investigated. Numerical simulations using the Wagner model were performed to explore the thickness distribution of the thermoformed products, and the results were compared with those of thermoformed products manufactured using industrial thermoforming machines.
The linear and nonlinear viscoelastic properties of the ABS polymers were studied by rheological measurements and the parameters of the constitutive model were obtained. In the SAOS experiments, the storage modulus and loss modulus were measured over a wide temperature and frequency ranges. A linear master curve was constructed to cover the wider range of frequencies. The discrete relaxation spectra and the parameters of the WLF equation were extracted from the master curve. To explain the nonlinear viscoelasticity, uniaxial extensional experiments were carried out and the parameters of the damping function were obtained. The extensional properties of the extruded ABS polymer sheets were found to significantly change depending on the orientation of the polymer chains. The type A specimen, which was oriented in the extrusion direction, had higher uniaxial extensional viscosity and stress than did type B, which was in the transverse direction. It was also found that, in type A, high resistance to material deformation and a wide range of uniform deformations are possible at 160°C and 170°C.
The nonlinear parameters α and β of the damping function were calculated by the optimization method. Using the Wagner model, numerical simulations were performed with the obtained nonlinear parameters, together with the viscoelastic properties. The thickness distribution of the numerical results was found to be qualitatively consistent with the products manufactured by an industrial thermoforming machine. Maximum difference between experiments and numerical results was observed at the edge in the upper compartment. Material deformations of the side walls were visualized, and grid deformations in the transverse direction in experiments were higher than those in the numerical simulations. Various reasons for the differences between the experimental and simulated predictions were analyzed with regard to the coefficient of friction, the orientation of the polymeric material, and the limitations of the constitutive model. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.