Modeling the Thermoforming Process of a Complex Geometry Based on a Thermo-Visco-Hyperelastic Model

: The thermoforming process is commonly used in industry for the manufacturing of lightweight, thin-walled products from a pre-extruded polymer sheet. Many simulations have been developed to simulate the process and optimize it with computer tools. The development of testing machines has simplified the simulation of this type of process, allowing researchers to characterize the behavior of the material at different temperatures and for large deformation to be closer to the real conditions of the process. This paper presents the results of a study on the modeling of the thermoforming process for an industrial demonstrator made from a high-impact polystyrene (HIPS) polymer. The HIPS shows a mechanical behavior that depends on the temperature and strain rate. In such conditions, a thermo-hyper-viscoelastic constitutive model is used to replicate the thermoforming process of the industrial demonstrator using ABAQUS/Explicit. Its behavior is determined via various experimental tests: uniaxial tensile tests at different temperatures and strain rates and Dynamic Mechanical Analysis (DMA). A comparison between the numerical and experimental results is carried out for the evolution of film thickness. The paper concludes with a discussion of possible improvements to be considered for future simulations of the thermoforming process using Abaqus, which presents complex challenges in terms of contact and material modeling.


Introduction
Thermoforming is a method of plastic parts manufacturing, in which a plastic sheet is heated and formed into a desired shape [1].The process is generally used in the packaging of consumer products and to manufacture large items such as automotive exterior parts or internal doors for refrigerators [2,3].The thermoforming process consists of four main steps: heating, pre-stretching, forming, and cooling (Figure 1).The sheet, or film, is heated by a radiative heater to a high enough temperature that it is malleable so it can be stretched onto a mold.The second step consists of creating a bubble by applying positive pressure.This step allows for decreasing the difference in the thickness between the edge and the center of the part.Then, the mold goes up until it touches the film.Vacuum pressure is used to remove the air from under the sheet.The part is cooled and ejected from the mold [4,5].
The forming of these materials is the subject of much research aiming to control the final properties of the products, such as thickness distribution.In this context, numerical simulation is generally used for the better control of these processes while avoiding expensive experimental tests.During the process, and due to the interaction between the polymer and the mold at different times, the thickness distribution is not uniform in the final part [6].The non-uniformity of such results can affect the part when the critical thickness occurs in a functional zone; therefore, manufacturers need to control and optimize their process to obtain the optimal distribution of the material [6].To simulate the rheological behavior of polymers for the thermoforming process, many studies have used various materials with various behaviors, like viscoplastic models, viscohyperelastic models, and visco-elastic [7][8][9][10][11].Takaffoli and al. [12] studied the influence of the temperature and the strain rate under the process conditions of a complex geometry.The authors used a thermal-mechanical finite element model to simulate the behavior of a polycarbonate sheet.Similar to the HIPS behavior, polycarbonate also depends on the temperature and strain rate; it also has a high stretchability above the glass transition temperature (Tg).To predict these characteristics, Takaffoli and al. developed a thermovisco-hyperelastic constitutive model based on the Arruda-Boyce model [13] and a powerlaw creep model.The thermal contact between the mold and the sheet was considered.The thickness distribution predicted by the adopted constitutive model was close to the experimental data.The investigation of O. Atmani et al. [14] was based on the simulation of the punch-assisted thermoforming of HIPS, for which they studied the influence of temperature and strain rate in thermoforming conditions.To predict the behavior of the material, a thermo-elasto-viscoplastic model derived from Gsell's behavior was used and implemented in Abaqus through a virtual user material subroutine, VUMAT.The model used takes into consideration the high impact of temperature, strain rate, heat exchange, and thermal contact.The results of punching force and sheet thickness estimated by the model were validated experimentally.A viscoelastic model was developed by Sánchez et al. [15], where the authors used the generalized Maxwell model.The model was implemented in Abaqus using Prony series parameters to replicate the thermoforming processes for polystyrene material under external loads and specific conditions of time and temperature.A good agreement between the numerical and experimental results was obtained.
Other studies have indicated that the thickness distribution of the final product can be affected by various factors, including the geometry of the clamping tools.Ertugrul et al. [16] emphasizedthe substantial impact of clamping ring geometry on wall thickness distribution.This underscores the importance of choosing clamping tool geometry following the product's geometry to attain a more uniform thickness distribution.Other articles explored the concept of size effects (SEs) in multiscale material processing and manufacturing, emphasizing the influence of changes in various influencing factors on the physical behaviors and phenomena of materials.In [17], Fu et al. provide an overview of size effects in micromechanical and microforming processes, highlighting the importance of understanding and managing these effects to ensure efficient manufacturing processes and high-quality products at the micro and meso scales.Recently, other authors have been interested in the application of digital twin technology in product design, manufacturing, and service, driven by big data.Tao et al. [18] present the methodology of digital twin-driven product design, manufacturing, and service, highlighting the need for high-performance computing, multi-physics/multi-scale interdisciplinary approaches, and deep learning.
The document also addresses the potential applications of digital twin technology in various phases of a product's lifecycle, including design, manufacturing, and service, with a focus on maintenance, failure analysis, and prediction.
In our study, a coupled thermal-mechanical finite element model is developed using Abaqus software.The thermo-visco-hyperelastic model, with a thermal contact conductivity at the interface of the mold and the sheet, is employed.This model takes into consideration the material's behavior depending on the temperature and strain rate.The latter is also implemented in Abaqus to simulate the thermoforming process of our demonstrator under industrial conditions.The paper particularly focuses on the experimental validation of the simulation on the final part's thickness distribution.

Thermo-Visco-Hyperelastic Material Model 2.1. Hyperelastic Model
Hyperelastic constitutive models are mathematical models that attempt to simulate the behavior of materials whose stress-strain relationships are nonlinear [19].These models are described by the deformation energy density function W, which is calculated as a function dependent on different variables associated with the deformation field and the material constants, as shown below: where λ 1 , λ 2 , λ 3 are principal stretches, I 1 , I 2 , I 3 are the invariants of a left Cauchy-Green strain tensor, and B and F are the strain gradient tensors [20].For hyperelastic models, the Cauchy stresses are calculated as follows: where J is the determinant of the strain gradient tensor F.
In Abaqus software, the strain energy potentials of incompressible isotropic elastomers can be approximately modeled with several forms: • Mooney-Rivlin model [20,21]: this model is highly recommended in previous existing studies for the simulation of material-forming processes.There exist many versions of this model based on the order of level, for instance: two-parameter, three-parameter, etc.The form of the strain energy for the two-parameter model is • Ogden model [20,21]: this model is also one of the most adopted hyperelastic models in the literature.The strain energy function is calculated as follows: where N, µ i , and α i are material parameters and D is a material constant related to the bulk modulus; • Marlow model [21]: this is a particular model, since it does not require any parameter calibration.Rather, a simple test is enough to ensure that the model follows the given data [22].The form of the Marlow strain energy potential is where W dev , W vol are, respectively, the deviatoric and volumetric parts.

Limitations of Hyperelasticity
Hyperelastic models are easy to use and calibrate in order to predict the behavior of rubber-like polymers.However, they cannot take into account the effect of strain rate [23].To capture strain rate effects, a viscoelastic model with hyperelasticity can be used.

Experimental and Numerical Tests
To study the effect of temperature and strain rate on the behavior of HIPS, several types of experimental tests were performed: DMA and uniaxial tensile tests.These tests were carried out below and above the Tg temperature.The choice of low-and hightemperature characterization tests was made to predict the behavior of the material at different temperatures.This allowed for simulating the real conditions during the experimental thermoforming process, taking into account the variation of the sheet temperature when it is in contact with the mold.Figure 2 shows the storage and loss moduli curves versus temperature.The loss modulus increases with the decrease of the storage modulus until reaching a maximum point, and then both curves decrease.The Tg temperature value is equal to 100 °C.This value is obtained using the peak of the loss modulus curve, as illustrated in Figure 2.

Tensile Tests
Tensile tests were conducted using the electromechanical testing machine of type INSTRON 5960 equipped with a 30 kN load cell.Displacement measurements were carried out using image correlation techniques (VIC 3D).The tensile experiments were performed at various temperatures and loading speeds.For each uniaxial tensile test, three dumbbellshaped specimens were selected.The dimensions of the specimen are specified in Figure 3.The diverse testing conditions and parameters employed are summarized in Table 1. Figure 4 shows true stress-true strain curves obtained at different temperatures; these curves illustrate the temperature effect on the HIPS behavior.For temperatures below Tg, the material has a glassy behavior and deforms until failure.For temperatures above Tg, the material presents a rubbery behavior.

The Constitutive Model Used in ABAQUS
To replicate the thermo-hyperelastic behavior of the HIPS material in ABAQUS, several models are available, including the Mooney-Rivlin, Ogden, and Marlow models.The selection between these models involved simulating the uniaxial tensile test on ABAQUS using different models.A correlation was established between the experimental and numerical tests to compare the stress-strain evolution.After testing several hyperelastic models, the Marlow model was selected due to its ability to accurately predict the material behavior at temperatures both below and above Tg (the glass transition temperature).The Marlow model assumes that the strain energy potential is independent of the second deviatoric invariant I 2 .This model is defined by providing test data that characterize the deviatoric behavior.In ABAQUS, a strain energy potential is constructed based on the provided test data to precisely reproduce the material's behavior, as shown in Figures 6  and 7.The input data used to reproduce the material's behavior as a function of temperature are shown in Table 2.
As previously mentioned, hyperelastic models do not consider the effect of the deformation rate.In this simulation, a hyper-viscoelastic coupling is used to implement the deformation rate effect.In Abaqus, viscoelasticity can be modeled for large-strain applications as a function of reduced time for time-dependent analyses [24,25].The viscoelastic model is defined in Abaqus using the Prony series, where the relaxation moduli are given as follows [26,27]: where G ∞ is the long-term modulus once the material is relaxed and τ i are the relaxation times.The parameters G ∞ , G i , and τ i have been determined through inverse identification using experimental data from DMA tests.Figures 8 and 9 demonstrate a good correlation with the experimental and numerical data for loss modulus and storage modulus.The Prony series parameters, identified as a result of this correlation, are provided in Table 3.These parameters will subsequently be utilized as input data for the FEM model.

Experimental-Numerical Modeling of the Industrial Demonstrator of the Thermoforming Process
The thermoforming process of the HIPS sheet of 1.4 mm thickness was simulated numerically using the thermo-visco-hyperelastic model.As shown in Figure 10, the numerical model is made of three parts: a mold, a HIPS sheet, and a support to hold the sheet.The mold and the support were modeled as isothermal rigid bodies, incapable of deformation but capable of contact and thermal interaction with the sheet, allowing conduction in the presence of contact and even a small gap [27].The HIPS sheet is meshed with shell elements (S4RT).The rectangular sheet was blocked on its outboard edges, as shown in Figure 11, and a very low pressure (0.1 bar) was applied on its surface, as shown in Figure 12, for the pre-stretching step.Then, a displacement was imposed on the mold and, finally, a negative pressure was applied to the sheet so that it took the shape of the mold.The Tables 4 and 5 provide a summary of all the boundary conditions and loads.

Contact Type Parts Properties
Mechanical contact Sheet and mold Hard contact and friction coefficient (f = 0.3) Thermal contact Sheet and mold Gap conductance 1 1 The data are given below.
The initial temperature of the sheet was defined using the experimental values obtained using a thermal camera during the experimental test, as shown in Figure 13.Considering the real conditions of the film temperature (Figure 14), a maximum temperature was applied on the red zone of 140 °C and an average minimum temperature of 100 °C was applied on the green zone.The temperature of the mold was also imposed; it was assumed that the temperature would remain constant (T mold = 50 °C).In addition, a mechanical contact was applied using a friction coefficient as well as a thermal contact using a gap conductance as mentioned in the Table 6.

Baseline Simulation
Figure 15 shows the different baseline simulations of the thermoforming process.Firstly, the HIPS sheet was formed into a bubble by applying uniform pressure (1-Prestretching) to improve the uniformity of thickness distribution.Then, the mold moved up until it attained its final position (two-mold moving).Finally, a negative pressure was applied on the surface of the HIPS sheet to push it against the mold (three-vacuum pressure) and to take its shape.

Temperature Distribution
Figure 16 shows the temperature distribution during the various stages of the simulation.During the pre-stretching step, the sheet is not in contact with the mold, and as a result, its temperature has not yet changed.It remains consistent with the initially imposed temperatures.For the second step, the mold comes into contact with the sheet.Due to the thermal contact, the temperature of the areas of the sheet in contact with the mold decreases.After the vacuum step is completed, all the surfaces of the sheet take on the shape of the mold, resulting in a decrease in temperature.

Stress Distribution
Figure 17 shows the von Mises stress distribution during the various stages of the simulation.The distribution of von Mises values shows a stress concentration in the first step due to boundary conditions with the fixed support.Then, for steps 2 and 3, there is a stress concentration in the areas that are most stretched by the rise of the mold.

Thickness Distribution
To validate the implemented numerical model, a comparison between the numerical and experimental results of the film thickness distribution was carried out.The thickness measurement points are presented below.Figure 18 shows the measurement points for the thickness of the study piece.A total of 125 pieces were used for measuring the thickness at various points.The measurement was carried out using a 3D imaging measurement machine, where a robotic measuring arm probed the piece to perform the measurement.Subsequently, an average measurement was considered for each point.Figure 19 shows the variation of the numerical thickness on one of the lateral faces (Figure 19) of the shaped part compared to the experimental data.For the measurement points 8, 9, 10, and 11, the average percentage error was around 13%.For the rest of the points, the average relative error was 30%. Figure 20 shows a comparison between the numerical and experimental results of the measurement points on the upper surfaces of the sheet, as shown in Figure 16c.From point 12 to 20, the average percentage error is about 14%. Figure 20 shows the evolution of the thickness for the rest of the points (50, 51, 52, 53, and 54) with an average relative error of 7%. Figure 21 illustrates the evolution of the thickness for the other lateral faces of the deformed shape.In this part, the difference between the numerical and the experimental results is the highest, for an average percentage of 48%.These results show that the error between the numerical and experimental results is lower (7%) at the upper surfaces of the deformed part.In this zone, the elements are less deformed and the temperature decreases faster compared to the evolution of the temperature in the lateral faces due to the heat transfer between the mold and the sheet.The elements at the lateral faces show a large deformation at a very high temperature, which makes the simulation more difficult and causes some stability problems.The results of our simulation are in agreement with those reported in the literature.In reference [12], a comparison was made between the numerical and experimental results for the thickness distribution of a PC sheet.A maximum error of 45% was mentioned.

Conclusions
In this paper, a finite element model was used to simulate the thermoforming process of an industrial demonstrator made of HIPS polymer.The HIPS behavior was evaluated by experimental uniaxial tensile tests and DMA tests.To implement the behavior of HIPS in ABAQUS finite element software, a visco-hyperelastic model, which allows for taking into account the effect of the strain rate and the temperature, was chosen.The Abaqus/Explicit solver was employed to simulate the process, where the thermal contact between the mold and the polymer sheet was included.To validate the proposed numerical model, a comparison was also carried out regarding the thickness distribution for several sections of the sheet.It was noticed that the visco-hyperelastic constitutive model can be simulated successfully.However, the nature of analysis contains both material and geometric nonlinearity, which could lead to considerable instability problems.The discrepancy between the numerical and the experimental thickness values is due to the following limitations of the constitutive model that we recommended are addressed in future works:

•
The type of mesh used in this study does not allow for using automatic selective mesh refinement to improve the deformation of a complex shape.

•
The influence of the friction coefficient on the evolution of the thickness was not studied.

•
Another important factor that was not studied is the conduction coefficient.

Figure 1 .
Figure 1.Thermoforming process: (a) heating: the sheet is heated using radiant heat until it reaches the optimal temperature for thermoforming; (b) pre-stretching: the sheet is mechanically stretched before molding, ensuring uniform thickness; (c) mold up: the mold goes up until it touches the film; (d) vacuum on: vacuum pressure is used to remove air from under the sheet and ensure close contact with the mold contours.
3.1.DMA Tests Dynamic Mechanical Analysis (DMA) is a technique used in materials science and engineering to analyze the mechanical and viscoelastic properties of materials as a function of temperature, time, and frequency.In this study, two types of DMA tests were conducted: • Temperature sweep: A frequency of 1 Hz was used, and a temperature ramp of 1 °C per minute was applied; • Frequency sweep: Frequencies ranging from 0.1 Hz to 10 Hz were used, with 10 points per decade.

Figure 3 .
Figure 3.The dimensions of a dumbbell-shaped specimen for the tensile test.

Table 1 .
Uniaxial tensile tests as a function of temperature and loading speed.

Figure 4 .
Figure 4. Experimental tensile curves of HIPS at different temperatures and a velocity of 0.15 mm/s.

Figure 5
Figure 5 illustrates the effect of strain rate for a temperature test of 60 °C.It should be noted that the stress level increases with the velocity of the load.

Figure 5 .
Figure 5. Experimental tensile curves of HIPS at different velocities and temperature tests of 60 °C.

Figure 6 .
Figure 6.Comparison of experimental and numerical tensile curves at different temperatures above Tg.

Figure 7 .
Figure 7.Comparison of experimental and numerical tensile curves at temperature test below Tg.

Figure 8 .
Figure 8. Correlation of experimental and numerical loss modulus curves.

Figure 9 .
Figure 9. Correlation of experimental and numerical storage modulus curves.

Figure 10 .
Figure 10.Different parts of assembly.

Figure 14 .
Figure 14.Distribution of experimental temperature values of the plate obtained through an IR camera.

Figure 15 .
Figure 15.Steps of the numerical simulation.(a): the sheet is mechanically stretched before molding, ensuring uniform thickness.(b): the mold goes up until it touches the film.(c): vacuum pressure is used to remove air from under the sheet and ensure close contact with the mold contours.

Figure 16 .
Figure 16.(a) Distribution of sheet temperature during the pre-stretching step, (b) distribution of sheet temperature during the mold up step, and (c) distribution of sheet temperature during the vacuum on step.

Figure 17 .
Figure 17.(a) Distribution of von Mises stress during the pre-stretching step, (b) distribution of von Mises stress during the mold-up step, and (c) distribution of von Mises stress during the vacuum on the step.

Figure 19 .
Figure 19.(a) Thickness distribution on side 1.(b) Experimental measurement points of the thickness on side 1.

Figure 20 .
Figure 20.(a,b) Thickness distribution on side 2. (c) Experimental measurement points of the thickness on side 2.

Figure 21 .
Figure 21.(a) Thickness distribution on side 3. (b) Experimental measurement points of the thickness on side 3.

Table 2 .
Material input data as a function of temperature used in FEM.

Table 3 .
Prony series parameters for FEM model.