Experimental investigation and constitutive modelling of the deformation behaviour of high impact polystyrene for plug-assisted thermoforming

This paper concerns the experimental and numerical study of the plug-assisted thermoforming process of high impact polystyrene (HIPS). The thermomechanical properties of this polymer were characterized at different temperatures and deformation rates. To study the influence of different parameters in the real conditions of plug-assisted thermoforming process, we carried out “plug-only” tests at different temperatures and plug velocities. To model the deformation behaviour of HIPS, we proposed a thermo-elastic-viscoplastic model, which we have implemented in Abaqus software. A thermo-dependent friction model was also proposed and implemented in Abaqus software. The parameters of the proposed models were identified by the inverse analysis method in the real conditions of plug-assisted thermoforming. The proposed models were validated with “plug-only” tests and plug-assisted thermoforming of yogurt container.


Introduction
Thermoforming is the process of heating a thermoplastic sheet to its softening point, pressing into or stretching it over the mould using vacuum pressure, air pressure, or mechanical force, and holding it in place while it cools and solidifies into the desired shape. Thermoforming is commonly used for food packaging, plastic toys, cafeteria trays, shower enclosures, and electronic equipment. A variety of thermoplastic materials can be used in this process, including Low Density Polyethylene (LDPE), High Density Polyethylene (HDPE), Polypropylene (PP) and High Impact Polystyrene (HIPS). HIPS is a lightweight and inexpensive thermoplastic often used in thermoformed packaging such as food and pharmaceutical containers [1]. In the plug-assisted thermoforming, a plug is used to force the softened material into the mould to facilitate the stretching of the sheet and to bring enough material in the bottom of the final product [1,2]. In most cases, the plug will not push the material completely into the mould but only part way. A pressure is then applied to draw the material against the cavity walls and complete the thermoforming operation. The main advantage of plugassisted thermoforming is its ability to give a better wall thickness uniformity than can be obtained without plugassist, especially for deep shapes. Several studies showed that the plug phase is the most important step in plugassisted thermoforming [2][3][4][5]. The final plug displacement corresponds to 70%-95% of the depth of the final product.
Several parameters control the quality of the final product: sheet material, plug material, sheet and plug temperatures and plug velocity [6,7]. The number of parameters and the complexity of the process require costly experimental studies [8]. To reduce experimental tests, numerical simulation is increasingly used to optimize the thermoforming process [2,9]. Several studies have shown that in plug-assisted thermoforming process, three main phenomena govern the quality of the final product in a coupled way: the material behaviour, sheet-plug contact and thermal transfers [1,2,10,11].
To model the behaviour of polymers in the plug-assisted thermoforming process, hyperelastic, viscoelastic and viscoplastic models are generally used. Dong et al. [9] used the Mooney-Rivlin hyperelastic model [12,13] to simulate the air pressure thermoforming of acrylic polymeric sheet, where the "perfect stick" condition is assumed between the mould and the polymeric sheet to obtain satisfying simulation results. Bernard et al. [1] used three hyperelastic models: Mooney-Rivlin [12,13], Arruda-Boyce [14], and Yeoh [15] to model the rheological behaviour of polystyrene in the plug-assisted thermoforming process. They showed that the hyperelastic approaches do not give satisfactory results and suggest using temperature and strain rate dependent elastoviscoplastic models. O'Connor et al. [16] developed a large strain thermally coupled (LSTC) material model for plug-assisted thermoforming of polypropylene yogurt container. They confirmed the requirement for the model to be fully thermomechanically coupled, and the process is most sensitive to plug-sheet contact friction and heat transfer parameters. Kittikanjnaruk and Patcharaphun [17] used commercial simulation package (T-SIM) to investigate the wall thickness distribution of HIPS and amorphous polyethylene terephthalate (A-PET) thermoformed parts. A nonlinear time-dependent K-BKZ model which was proposed by Kaye [18] and Bernstein et al. [19] was used. Hosseini et al. [20] proposed a nonlinear viscoelastic model with new strain energy function to model the plug-assisted thermoforming of ABS sheets. Most of these models are valid on a limited range of temperatures and strain rates or developed for a specific material. Wang et al. [21] developed a finite-strain and thermomechanically-coupled constitutive model, which is based on a tripartite decomposition of the deformation gradient into elastic, viscoplastic, and thermal components for the simulation of gas-blow forming of amorphous glassy polymers.
Besides the thermomechanical behaviour of the material, sheet-plug interactions have a very important effect on thickness distribution of the final product [1,3,20]. Nam et al. [22] investigated the effects of mould wall surface by examining a perfect slip and no slip conditions. They found a clear difference between the two conditions and concluded that further work is required in the area of heat transfer between the mould and the sheet. Petty [23] stated that a friction model should ideally relate to the instantaneous conditions that exist at all contact points in the model. By studying the effects of contact conditions in plug assisted thermoforming, Collins et al. [1] concluded that contact conditions (friction and heat transfer) must be included in the simulations. Martin et al. [24] carried out experimental measurements of the polystyrene and polypropylene friction with two different types of plug. The authors showed that sheet-plug friction is dependent on the temperature and the material of the plug. Morales et al. [11] measured the friction coefficient of HIPS under different conditions. They showed that the friction coefficient increased with the sheet temperature, but the variations were limited regarding the plug material or the plug velocity. Laroche et al. [25] fabricated rounded containers from PP sheet using Delrin (Acetal and POM) and syntactic foam plugs. Friction coefficient between the polymer sheet and the plug was measured as a function of temperature and part thickness distributions were compared with predictions. The authors concluded that the surface roughness, the temperature of the interface, the contact pressure and the velocity of the slip also affect the contact friction.
In view of this literature review, to properly perform simulation of the plug-assisted thermoforming process, one needs a thermomechanical material behaviour model covering stain rates and temperatures ranges of thermoforming, an accurate modelling of the sheet-plug contact tacking into account the heat transfers and the plug material, and relevant experimental tests to calibrate the proposed model parameters and their coupling. A lack of research concerning the deformation in more than one axis is evident. To address these issues, firstly, we conducted experimental study based on different mechanical tests to represent the deformations encountered during plugassisted thermoforming. The thermal dependence of the elastic behaviour was analysed using dynamic mechanical analysis (DMA). Uniaxial tensile tests at different temperatures and strain rates provided a broad experimental description of intrinsic behaviour of the material. The plugonly tests were used to give a more realistic behaviour of the complex deformation state encountered under and around the plug and during the plug phase. Secondly, we proposed a phenomenological modelling of the HIPS behaviour based on a thermoelastoviscoplastic formalism which was implemented in the finite element simulation code Abaqus via VUMAT subroutine. The proposed model was inspired from G'Sell et al. model [26][27][28][29] which is a phenomenological elastoviscoplastic law widely used to model the mechanical behaviour of different polymers [30][31][32][33][34]. The effect of the contact between the sheet and the plug was considered by the temperature dependence introduced into the Abaqus simulation code via VFRICTION subroutine. The parameters of the proposed models are then identified using an inverse identification method. The proposed methodology was finally validated by comparing experimental and finite element simulation results of plugassisted thermoforming of a dairy packaging.

Material and experimental testing 2.1 Material
In this study, we have used commercial extruded high impact polystyrene (HIPS) sheets with 1.1 mm thickness supplied by Intraplás Ind ustria Transformadora de Plásticos S.A. (Portugal). HIPS has a higher impact strength than homopolymer PS. This type of PS is produced by adding around 5%-10% rubber or butadiene copolymer. This increases the toughness and impact strength of the polymer and results in a very stiff product used in packaging applications for a wide range of applications in the food, medical, consumer goods, cosmetics, industrial and horticultural markets among others. HIPS is an amorphous polymer which is a translucent white in its natural state. This polymer is compatible with food packaging applications and has inherent moisture barrier properties that make it suitable for short shelf life products.

Tensile tests
Uniaxial tensile tests were carried out at different temperatures and crosshead velocities using Instron (33R 4204) testing machine equipped with 1 kN loading cell and an eclectic oven. For this purpose, samples have been prepared from the HIPS sheets according to ISO 527 standard test (Type 1 sample).
The stress-strain typical curves obtained from tensile tests at different temperatures and crosshead velocities are shown in Figure 1. They exhibit typical behaviour of glassy polymers [35][36][37] where three stages can be distinguished. At low strains a viscoelastic response can be observed, followed by yielding behaviour and finally plastic behaviour at high strains. The viscoelastic response can be divided in a linear viscoelastic part at low stress and a nonlinear viscoelastic region at higher stresses.
At a crosshead velocity of 500 mm/min a significant change in the curves is observed between the different temperatures, especially between 100°C (Fig. 1b) and 110°CÀ120°C (Fig. 1a). The temperature range [100°C-120°C] is characterized by the transition from a ductile glassy behaviour for T 100°C to a viscous behaviour for T ≥ 120°C. This behaviour can be explained based on the molecular motion of polymer chains. At low temperature, the polymer structure is rigid and compact. Therefore, the molecular motion in HIPS is small because of the low kinetic energy of molecules, leading to high elastic modulus and tensile strength. As the temperature increases, the kinetic energy of molecules increases, which increases the molecular motion resulting in a decrease in the mechanical properties. Near the glass transition temperature (100°C-110°C), the molecular motion is further enhanced for viscous flow and the samples become soft and rubbery, leading to very low values of elastic modulus and tensile strength.
From the results shown in Figure 1, we can extract both modulus and yield stress at different experimental conditions. Using the Ree-Eyring representation [38], the modulus and yield stress are plotted as a function of logarithm of strain rate at different temperatures in Figure 2a, showing a linear dependence versus the decimal logarithm of the strain rate for both properties. We have also plotted the modulus and yield stress as a function of temperature at different strain rates in Figure 2b.
There is a significant change in the slopes of the plots after 100°C corresponding to the a transition zone.

DMA tests
Dynamic Mechanical Analysis (DMA) is a technique used to measure the mechanical properties of a wide range of materials. Dynamic Mechanical Analysis test was carried out in a DMA (Q800, TA Instruments) using dual cantilever bending mode. Rectangular shaped samples with 55 mm length and 15 mm width were cut out from the sheets using a punch. These samples were tested in a temperature range from 20°C to 140°C at a heating rate of 3°C/min and a frequency of 1 Hz. In dual cantilever bending mode, the sample is clamped at both ends and flexed in the middle. Figure 3 shows a scan of HIPS plotted as the storage modulus (E 0 ), loss modulus (E 00 ) and tan d versus temperature. The strong drop in storage modulus observed above 95°C is related to the a relaxation. A common measurement on polymers is the glass transition temperature T g that can be measured by the storage modulus onset point, by the loss modulus peak, or the peak of tan d. The loss modulus peak is considered in this study leading to a value of T g ≈ 105°C.
Furthermore, polystyrene has b relaxation below the primary relaxation a around 76°C as shown in Figure 3. This b process was attributed to the local oscillation mode of backbone chains by Yano and Wada [39]. An NMR narrowing at around 77°C observed by Odaiba et al. [40] was also attributed to the b relaxation. The b peak is separated from the a peak only at low frequencies and merges into the a peak at high frequencies [39].

Plug-only tests
In this study, plug-only test at high velocities is used to study the mechanical behaviour of HIPS under multiaxial and large deformation loading allowing to be close to the process conditions. The plug-only test consists of hitting at  Figure 4. Influence of sheet temperature, plug velocity and plug lubrication was studied to reproduce HIPS behaviour under thermoforming conditions. The plug force as a function of the plug stroke curves obtained at different temperatures with an aluminium plug heated to the same temperature as the HIPS sheet and with plug velocity V = 1 m/s are shown in Figure 5. This figure shows a rigid-fragile behaviour of the HIPS sheet at temperature T = 110°C with low deformation before it breaks at 30 mm of plug stroke. At T = 120°C, the deformation of the sheet is greater, but its rigidity remains significant corresponding to a ductile behaviour. At highest temperatures (T = 130°C and T = 140°C), we see a rubbery behaviour in which the material undergoes significant deformations without rupture. The maximum plug force decreases significantly, and the slope of the curves become very small when the temperature increases.    Figure 6 shows the effects of plug velocity and plug lubrication on the plug force obtained at T = 130°C with a hot aluminium plug. An increase in the plug velocity results in an increase in the rigidity of the HIPS sheet and the force required to deform it (Fig. 6a). Figure 6b shows that a higher plug force is applied to the sheet without lubrication (130 N) compared to the use of the lubricant (100 N). The photos show that the thickness of the sheet in contact with the plug is lower when lubrication is used, whereas when the plug is not lubricated more material is present under the plug due to friction.
3 Modelling of the thermomechanical behaviour

Elastoviscoplastic constitutive model
In the present work, the elastoviscoplastic model proposed by G'Sell et al. [26][27][28][29] is adapted to describe the thermomechanical behaviour of HIPS by introducing the temperature dependence as follows: where e is the total strain, e p is the equivalent plastic strain, _ e is the strain-rate, E(T) is the elastic modulus, s 0 (T) is the yield stress, K(T) is a scaling factor, 1 À exp Àwe p À Á À Á is a viscoelastic term, exp hðT Þe p À Á takes into account the very important strain-hardening observed at large deformations. Finally, the term _ e m expresses the strain-rate sensitivity as a power law function.
Furthermore, the strain softening behaviour corresponding to the post-yield stress drop when the material initially enters the viscoplastic regime was addressed in G'Sell original model by the multiplicative term (1 + a . exp(À be)). This term was not considered in this study, because the thermoforming of HIPS is carried out at high temperatures where the strain softening is less important as shown in Figure 1. In addition, this term uses two more parameters to identify.
A sensitivity study shows that the two material parameters w and m are temperature-independent, while the parameters K(T) and h(T) depend on the temperature in the same way as the modulus E(T) and the yield stress s 0 (T) following an Arrhenius evolution:  where E 0 , A E , s 00 , A s , K 0 , A K , h 0 , and A h are material parameters determined from experimental results by the procedure indicated in the following sections. The elastoviscoplastic constitutive law was implemented in the commercial nonlinear finite element software Abaqus/Explicit through a VUMAT Fortran subroutine. Abaqus/Explicit provides the user subroutine VUMAT the strain increment for the current time step, the stress tensor and the temperature at the beginning of the current increment, the time increment, and a table of solution dependent state variables. The VUMAT subroutine computes and returns the value of the stress tensor and the state variables at the end of the increment for each integration point.

Friction model
In plug-assisted thermoforming, the temperature at the sheet-plug contact zone changes during the process. This temperature depends on many parameters: the polymer and plug temperatures, the plug velocity, and the polymer and plug materials. Several studies [1,3,5] showed that friction is highly related to this interface temperature. Laroche et al. [25] proposed the following exponential-type function to represent friction temperature dependence: fðT where f 0 , f 1 and f 2 are the coefficients of the friction model which can be identified from experimental results by the procedure indicated in the following sections. The friction model was implemented in the commercial nonlinear finite element software Abaqus/Explicit through a VFRICTION Fortran subroutine. Abaqus/Explicit provides the user subroutine VFRICTION the contact points and the temperature at the beginning of the current increment, the time increment, and a table of solution dependent state variables. The VFRICTION subroutine computes and returns the tangential frictional force at each contact point in local coordinates and the state variables at the end of the increment.

Simulation of plug-only test
The models implemented in ABAQUS/Explicit were applied to the simulation of plug-only test, where the material undergoes multiaxial deformation with high deformation rates. The rectangular sheet of dimensions 200 Â 200 Â 1.1 mm 3 is meshed with 5128 S3RT elements. The sheet is clamped between the mould and the countermould leaving free a circular central part of diameter 90 mm which will undergo the deformation as shown in Figure 7. Its temperature was set at the temperature of the simulated test. The rigid tools (mould and counter-mould) are discretized with isothermal rigid shell elements R3D4 and their temperature was set to 25°C. The plug was meshed with 322087 C3D4T 4-node thermally coupled tetrahedron, linear displacement and temperature. Its velocity corresponds to that of the simulated test.
The lower face of the plug is positioned 20 mm above the sheet. The plug was constrained as rigid body, but not as isothermal to consider heat exchanges with the polymer sheet. The conductive heat transfer between the contact surfaces is assumed to be defined by: where q is the heat flux per unit area crossing the interface from polymer sheet to plug surface, T Polym and T Plug are the temperatures on the polymer and plug surfaces respectively, and k is the gap conductance. The heat exchange of the polymer sheet with the external environment is also considered: where q is the heat flux per unit area, T Polym is the temperature of the polymer, T amb is the ambient temperature, and h is the global exchange coefficient (convection + radiation).

Identification procedure of model parameters
The procedure of identification for the constants of the elastoviscoplastic and friction models is summarized in the schematic diagram given in Figure 8, and can be described as follows: -Firstly, the temperature dependent elastic parameters were determined. Elastic moduli were determined from tensile and DMA tests, and yield stresses were extracted from tensile tests. The temperature dependent parameters used in equation (2) were then determined by fitting the experimental curves. -Then, a numerical inverse identification method is used to determine the viscoplastic model parameters.
The experimental plug force curves of a set of lubricated plug tests performed under specific conditions (temperature and plug velocity) were compared to the calculated curves obtained by the direct simulation of these tests. -Finally, a numerical inverse identification method is used to determine the friction model parameters. The experimental plug force curves of a set of unlubricated plug tests performed under specific conditions (temperature and plug velocity) were compared to the calculated curves obtained by the direct simulation of these tests by using the friction model.

Identification of elastic properties from DMA and uniaxial tensile tests
Calibration of the elastic properties was conducted using previously detailed DMA and uniaxial tensile tests. Firstly, elastic moduli were determined from tensile and DMA tests, and yield stresses were extracted from tensile tests. Then, the Arrhenius-type dependency of elastic modulus and yield stress as a function of temperature equation (2) were fitted with these experimental results.

Inverse identification of viscoplastic properties from plug-only tests
To determine the unknown set of viscoplastic parameters P = (K 0 , A K , w, h 0 , A h , m) in the constitutive equations (1) and (2), a numerical inverse optimization method is used. Considering a set of lubricated plug-only tests performed under specific conditions (temperature and plug velocity), the calculated plug forces depending on the unknown parameters P, were obtained by the simulation of these tests, and were compared to the measured results through the cost function defined as: where F sim (P, T j , V j , t i ) and F exp (T j , V j , t i ) are the simulated and measured plug forces at the t i sampling point on the jth test curve at temperature T j and plug velocity V j . M j is the total number of experimental curves, N i is the number of sampling points used for comparison for each experiment, and p j is a weighting factor (0 < p j 1). The minimization process of the cost function was driven by a NSGA-II algorithm (Non-dominated Sorting Genetic Algorithm II) [41]. The ABAQUS/Explicit code with the developed VUMAT subroutine was used to perform the plug-only test simulations and a python program was developed to extract the plug force curve. An in-house developed interface was used to transfer the successively updated sets of parameters P from NSGA-II algorithm to ABAQUS, and to calculate the values of the cost function according to equation (6) and transfer them to the NSGA-II algorithm.
The optimisation process begins by selecting an initial design population randomly distributed in the space design and then running the NSGA-II algorithm with 20 generations to keep the simulations computationally tractable. The best individual of the last generation is then used as a start point for the gradient-based Levenberg-Marquardt algorithm [42]. This iterative process runs until the convergence conditions are achieved.
To make provision against non-uniqueness of the solution, some measures have been taken. Firstly, the weighting factors 1/M j are introduced in equation (6) to ensure that the number of selected sampling points for each experiment doesn't affect the solutions. Secondly, the experimental data base covers multiple testing conditions: tensile tests at temperatures ranging from 80°C to 140°C and three crosshead velocities (50, 100 and 500 mm/min), plug-only tests at temperatures ranging from 110°C to 140°C and three plug velocities (0.2, 0.5 and 1 m/s) with and without lubrication, are taken into account to calculate the cost function given in equation (6). However, it is seen that some sets produce indistinguishable plug forces, indicating that it is difficult to eliminate the nonuniqueness completely when inversely identifying the material parameters from real experiments. Shrot and Bäker [43] came to the same conclusion during the inverse identification of Johnson-Cook material model for the simulation of the machining process.

Inverse identification of friction parameters from plug-only tests
In this step, the elastic and viscoplastic parameters were fixed, and the unknown set of friction parameters P = (f 0 , f 1 , f 2 ) in equation (3) were determined by using a numerical inverse optimization method. The calculated plug forces depending on the unknown parameters P, were obtained by the simulation of unlubricated plug-only tests using the friction model and were compared to the measured results through the cost function defined in equation (6). The same optimisation process as for identification of viscoplastic parameters is used.  The experimental results were fitted by the Arrhenius-type functions defined in equation (2) in the temperature range (100°C T 140°C) with good correlation coefficients (R 2 ≈ 0.93). For temperatures lower than 100°C, elastic modulus and yield stress are almost constant. The fitted parameters (E 0 , A E , s 00 , A s ) are given in Table 1.

Plug-only test results
Firstly, the simulations of the pushing tests were carried out by considering only the rheological behaviour of HIPS and eliminating the friction and the sheet/plug thermal transfer based on the experimental results of the plug-only tests with lubrication. Thus, the frictionless contact condition was used allowing the elimination of the friction between the polymer sheet and the plug, and the same temperature was considered for the sheet and the plug. The identified viscoplastic parameters P = (K 0 , A K , w, h 0 , A h , m) are given in Table 1. The experimental and predicted curves of plug force vs plug stroke without friction at different temperatures and plug velocities are shown in Figure 10. A good agreement between predictions (solid lines) and measurements is observed for the tested temperatures and plug velocities. However, the calculated forces at low temperatures and low velocities are underestimated compared with the measured results. Indeed, we have affected greater weighting factors (p j ) for high temperatures and high velocities in the identification process to favour the thermoforming conditions of HIPS.
For the identification of friction parameters, the simulations of the pushing tests were carried out by considering the rheological behaviour of HIPS with thermomechanical friction between sheet and plug based on the experimental results of the plug-only tests without lubrication. The identified friction parameters P = (f 0 , f 1 , f 2 ) are given in Table 2. The experimental and predicted curves of plug force vs plug stroke with friction at different temperatures and plug velocities are shown in Figure 11. A good agreement between predictions (solid lines) and measurements is observed for the tested temperatures and plug velocities. We have used the same approach as for frictionless identification which explains why the calculated forces at low temperatures and low velocities are underestimated compared with the measured results.
The comparisons of the curves of the plug force vs the plug stroke allow to verify the overall behaviour of the proposed models, but do not provide information about the local behaviour of the material. Thus, we have analysed the thickness distribution of the sheets at the end of plug-only process. Figure 12 shows photographs after plugging stage and the corresponding calculated thickness distribution obtained with and without friction at sheet   temperature of 130°C and plug velocity of 1 m/s. The experimental and simulation results show that without friction the material slides along the plug and the deformed sheet is thinner than with friction. This is confirmed by the comparisons of the experimental and numerical thickness profiles under different experimental conditions given in Figure 13. The thickness profile is reasonably well reproduced by the numerical prediction with parameters identified by overall identification procedure (solid lines).

Thermoforming simulation
In order to test the prediction capability of the numerically identified model, the experiment and FE simulation of the thermoforming of a yogurt container was performed. Figure 14 shows the experimental device for the plugassisted thermoforming used in this work. The mould is made of transparent polycarbonate on the sides, allowing to visualise the forming of the container, and of aluminium on its base allowing a rapid cooling of the container.
The polymer sheet temperature was set at 140°C and the plug and mould temperatures were set at ambient temperature. Same meshing as for the plug-only test was used for the sheet, the plug and the counter-mould (Fig. 7).
The mould was discretized with 7736 linear triangular isothermal rigid shell elements of type R3D3 (Fig. 15). The lower face of the plug is positioned 20 mm above the sheet and is subjected to an 80 mm stroke at a velocity of 1 m/s. At the end of the plug stroke, air is blown at a pressure of 5 bars to press the polymer sheet onto the mould. The total depth of the yogurt container is 67.5 mm and the HIPS sheet is deformed about 89% during the plugging step and about 11% during the blowing step. Figure 16 shows the thickness distribution at different stages of plugging (50%, 75% and 100%) and blowing (10%, 20% and 100%). We can notice that at the end of plugging, the minimum thickness is observed on the sides and the bottom of the container is thicker. The thickness profile along the centreline of the yogurt container obtained experimentally and by simulation are compared in Figure 17 showing good agreement between the two results.

Conclusions
In this work, we proposed a thermoelastoviscoplastic constitutive law and a friction model for the simulation of plug-assisted thermoforming process. The models were implemented in the ABAQUS/Explicit software using the VUMAT (constitutive law) and VFRICTION (friction model) user subroutines. For the identification of the model parameters a three steps methodology was proposed. The temperature dependent elastic material parameters were determined from experimental uniaxial traction and DMA tests. The temperature dependent viscoplastic parameters were determined using inverse optimisation procedure by comparing experimental and calculated plug force curves of a set of lubricated plug-only tests performed under specific temperature and plug velocity conditions. The friction model parameters were also determined by inverse identification method by comparing experimental and calculated plug force curves of a set of lubricated plugonly tests performed under specific temperature and plug velocity conditions. It was found that the sheet   temperature, the plug velocity and lubrication play significant role in governing the thickness distribution of the deformed sheet. Plug-only experiments validated the preponderant effect of plugging stage on the final thickness distribution. The experimental and simulation results show that without friction the material slides along the plug and the deformed sheet is thinner than with friction. Plugassisted thermoforming simulations were performed using the identified thermoelastoviscoplastic constitutive model for the polymer sheet and using temperature dependent model of friction between the plug and the sheet. The simulated thickness distribution showed qualitative agreement with experimental results.