# Finite Element Analysis of PMMA Stretch Blow Molding

Authors: Afef Bougharriou; Mohieddine Jeridi; Mohamed Hdiji; Anoir Boughrira; Kacem Saï
Journal: International Journal of Manufacturing Engineering
 (2014)

---

## Abstract

The electric bubbles are a useful product made of PMMA material. They are produced by the stretch blow molding process. Thickness, which reflects the quality of the electric bubble, is a crucial parameter that deserves special attention for the molding process. In this work, finite element simulations of the stretch blow molding process are performed aiming at the determination of the preform geometry to ensure homogeneous thickness of the finished product. The geometrical parameters of the preform are optimized allowing a better homogeneity thickness compared to existing data. The predicted parameters allow the improvement of the thickness distribution. The standard deviation of the thickness is reduced to about 95% compared to the existing bubble.


---

## Body

## 1. Introduction

Stretch blow molding is one of the most important manufacturing processes of polymer via injection process [1–3] or extrusion process [4, 5] to produce hollow plastic parts. The majority of clear bottles are produced using injection stretch blow molding process. Accordingly, several works were devoted to the simulation of the stretch blow molding in order to optimise the process and improve the quality of these containers. Chung [6] proposed a finite element (FE) model for injection stretch blow molding of PET bottles using the software ABAQUS. These simulations focused on the optimisation of the side wall thickness of the bottle by adjusting the variations of histories of plunger movement and gas pressure. Cosson et al. [7] used the so-called constrained natural elements method to simulate and optimize the stretch blowing of a PET preform by assuming an axisymmetrical and viscoplastic preform model under an internal pressure. The thickness of the bottle is linked to the variation of the stretch rod speed and time where the blow pressure is applied. Cosson et al. [8] simulated the stretch blow molding process of PET bottles at the usual process temperature. The simulation results modeled the crystalline microstructure evolution and predicted the elastic behaviour of the final bottle. It also provides Young’s modulus distributions in the bottles. Tan et al. [9] developed a 2D isothermal FE simulation of the injection stretch blow molding process for PET containers. To inflate the PET preform, two approaches were used in the simulation: (i) a direct pressure input and (ii) a constant mass flow rate of air. The results show that the simulation with the fluid flow method gives good prediction of the volume versus time curve and the preform shape evolution. Tan et al. [9] concluded that applying the pressure via a mass flow rate of air is the most appropriate method to model the pressure inside the preform. Salomeia et al. [10] simulated the stretch blow molding process for PET bottle by the means of thermodynamic model implemented in an ISBM. The simulations provide the evolution of the pressure and the temperature profiles on both inside and outside the surfaces. Yang et al. [11] proposed a FE model for the injection stretch blow molding process of PET bottles aiming at optimizing the process conditions (i.e., the bottle wall thickness, the distributions of the material, the temperature, and the strain). Gupta et al. [12] simulated the blow molding process for PET bottle in order to study the thickness and the stress distribution. Zimmer and Stommel [13] presented FE simulations for stretch blow molding with free blow trials of a PET bottle. The model provides the temperature distribution and the stress-strain fields. The poly(methyl methacrylate) (PMMA), material of this study, is an economical, versatile general purpose material. Various types of PMMA presenting high mechanical properties are used in different fields and applications. Numerous studies were devoted to the modeling of the mechanical behavior of this material; see, for instance, [14–17]. The FE studies related to the simulations of stretch blow molding process for PMMA material are scarce. Dong et al. [17] presented a FE simulation of thermoforming process of rectangular bubble inflation. The obtained critical material parameters of PMMA have been used in the simulation of free inflation of a bubble profile. The light bubbles are highly demanded in the Tunisian market. The Tunisian market imports these bubbles or the manufactured molds from abroad. IMR, a Tunisian company specialized in design and manufacture of injection tools, aims at producing these bubbles at lower cost and good quality. IMR needs to optimize the preform geometry leading to quasihomogeneous distribution of the thickness along the bubbles. In this work, FE simulations of the stretch blow molding process are performed aimed at the determination of the preform geometry to ensure homogeneous thickness of the finished product. The paper is organized in the following manner: in Section 2 the preform geometry (i.e., mesh), the boundaries conditions, the loading, and the material behavior are detailed. Section 3 is devoted to the discussion of the FE results in terms of stress and bubble thickness distributions. Finally, recommendations are given on the geometric design of the preform.


## 2. Finite Element Modeling

The simulations are performed using the FE software ABAQUS. The material behavior may affect in a considerable way the obtained results. To this end, the FE model of the molding process is investigated taking into account the hyperelastic behavior of the PMMA at 150°C by the means of Mooney-Rivlin model.


### 2.1. Geometry and Meshing

The design parameters of the preform are shown in Figure1(a). The geometry is controlled by the three independent parameters: e, α, and H. Assume that the deformation of the preform occuring under constant volume (V
=
V
0) reduces the number of parameters from 3 to 2 as will be shown below. The FE analysis is performed using six-noded axisymmetrical elements with a full integration scheme. The preform is designed according to the pairs (
α
,
e
), Table 1. The mesh size is chosen in such a way that subsequent refinement does not affect the stress-strain behavior much. A typical FE mesh is shown in Figure 1(b).Table 1
Number of elements in the FE model.
α
/
e
3.6
3.8
4
4.2
2.1
4221/1872
4005/1776
3825/1696
3609/1600
2.3
—
4401/1952
4149/1840
3879/1720
2.5
—
—
—
4329/1920Geometry of the preform.
(a)
Geometrical parameterization
(b)
Typical mesh (zoom)The volume of the preform is calculated as the sum of a semisphere and a conical cylinder as shown in Figure1. After some algebraic manipulations, the expression of the volume of the preform is given by the three terms A, B, and S:
(1)
V
=
A
+
S
-
B
with
(2)
A
=
π
3
H
′
(
R
1
2
+
R
2
+
R
1
R
)
,
B
=
π
3
H
′
(
R
2
2
+
r
2
+
R
2
r
)
,
S
=
4
π
6
(
R
3
-
r
3
)
.
A, B, and S depend on the geometrical characteristics of the preform:
(3)
R
=
D
-
(
H
′
-
K
)
tan
α
,
r
=
(
D
-
e
)
-
(
H
′
-
K
)
tan
2
,
R
1
=
(
D
+
K
tan
α
)
,
R
2
=
(
D
-
e
+
K
tan
2
)
,
H
=
R
+
H
′
.
It is assumed that the deformation of the preform occurs under constant volume during the mold flow process (V
=
V
0). This leads to a third degree equation in e:
(4)
V
=
V
0
⟹
e
3
+
a
e
2
+
b
e
+
c
=
0
,
where
(5)
a
=
-
3
(
D
-
(
H
′
-
K
)
tan
2
)
-
3
2
H
′
,
b
=
3
(
D
-
(
H
′
-
K
)
tan
2
)
2
+
3
2
H
′
(
2
D
+
(
2
K
-
H
′
)
tan
2
)
,
c
=
3
2
π
(
A
-
V
0
)
+
1
π
(
R
3
-
(
D
-
(
H
′
-
K
)
tan
2
)
3
)
-
1
2
(
(
D
+
K
tan
2
)
2
+
(
D
-
(
H
′
-
K
)
tan
2
)
2
1
2
+
(
D
+
K
tan
2
)
(
D
-
(
H
′
-
K
)
tan
2
)
1
2
)
.
Analytical solution can easily be obtained for the third degree (4). Two cases are distinguished depending on the value of the determinant Δ
=
q
2
+
4
p
3
/
27 (q
=
2
a
3
/
27
-
a
b
/
3
+
c and p
=
b
-
a
2
/
3). For the studied problem Δ is negative, so that 3 real zeros verify (4). Each solution k (k
=
0, k
=
1, and k
=
3) is given by the means of polar coordinates:
(6)
e
=
2
r
1
/
3
cos
⁡
(
θ
3
+
2
k
π
3
)
with r
=
0.5
q
2
-
Δ and θ
=
arctan
(
-
-
Δ
/
q
). Only one solution lies on the allowed range. The diagram in Figure 2 can be divided into regions. The divisions between these regions are bounded by the curve with the equation H
=
H
max
⁡. The feasible region is above this curve, whereas the infeasible region below the curve leads to preform geometry in which the height H is greater than the height of the mold.Figure 2
Limit curve.


### 2.2. Boundary Conditions and Loading

The following boundary conditions are considered for the FE analysis.(i)
The “container” is fully constrained and the preform could move only along the axial direction.
(ii)
A displacement along thez-axis (Figure 3) to the bottom of the mold is applied first to the preform. This displacement is preconized by the manufacturer. Accordingly, the preform moves down to be in contact with the mold.
(iii)
A pressure of 0.3 MPa was applied in the second step to the inside of the preform with a gradual increase.
(iv)
A friction law, corresponding to Coulomb law with a constant friction coefficientμ
=
0.1, is considered between the preform and the “container.” It is worth noting that some simulations (see Section 3) in which μ was taken in the range of 0.1 to 0.5 have shown that the results were not affected by the value of the friction coefficient.Figure 3
Loading.


### 2.3. Material Behavior Law

To describe the mechanical behavior of the PMMA, hyperelastic and viscoelastic models are the two commonly used approaches. The emphasis of current study is a phenomenological hyperelastic approach due to its better suitability for the PMMA which is isotropic, highly elastic, and incompressible for the forming conditions [17]. It is already established that the Mooney-Rivlin model is accurate in simulating polymer sheet forming; see, for instance, [17, 18]. It should be noted that this approach fails to represent the speed influence. The material behavior is described with a specific free energy density W which is a function of strain and temperature. For isotropic behavior, W depends on the invariants of the transformation tensor F which is given by
(7)
d
x
=
F
d
X
,
B
=
F
T
F
,
where d
X is the elementary vector in the initial configuration, d
x is the transformation of d
X in the actual configuration, F is the transformation gradient tensor, and B is the Cauchy-Green Gauche tensor. The stresses can be obtained as the partial derivatives of the strain energy density function with respect to invariants I
1, I
2, and I
3. For isotropic materials, W is a linear combination of the invariants I
1, I
2, and I
3 of the strain tensor W
=
W
(
I
1
,
I
2
,
I
3
), where I
1
=
trace
B, I
2
=
(
1
/
2
)
(
(
trace
B
2
)
+
(
trace
B
)
2
), and I
3
=
det
⁡
(
B
). These invariants can be expressed with the elongations (λ
i, i
=
1,2, and 3) as follows:
(8)
I
1
=
λ
1
2
+
λ
2
2
+
λ
3
2
,
I
2
=
λ
1
2
λ
2
2
+
λ
1
2
λ
3
2
+
λ
2
2
λ
3
2
,
I
3
=
λ
1
2
λ
2
2
λ
3
2
.
The principal stretch ratio λ is the ratio of the gauge lengths after and before the test (λ
=
L
/
L
0). The hyperelastic theory (Mooney-Rivlin model) is adopted in this study for the simulation of the nonlinear material behavior. The strain energy density W is expressed by Mooney-Rivlin model as follows:
(9)
W
=
C
10
(
I
1
-
3
)
+
C
20
(
I
1
-
3
)
,
where C
10 and C
20 are two material parameters. For the case of incompressible material under uniaxial elongation (λ
1
=
λ, λ
2
=
λ
3
=
1
/
λ), the true Cauchy stress differences can be calculated as
(10)
σ
11
-
σ
33
=
2
C
10
(
λ
2
-
1
λ
)
+
2
C
20
(
λ
-
1
λ
2
)
,
σ
22
-
σ
33
=
0
.
In the case of simple tension, σ
22
=
σ
33
=
0, the uniaxial behavior law writes
(11)
σ
11
=
(
2
C
10
+
2
C
20
λ
)
(
λ
2
-
1
λ
)
.
The Mooney-Rivlin parameters of PMMA material are taken from the literature. One element test subjected to monotonic tensile loading is carried out in order to check the accuracy of the material parameters in simulating the experimental stress/elongation curve (Figure 4). It can be seen that the experimental data are well described by the Mooney-Rivlin model for the parameters values of C
10
=
-
0.0303 MPa and C
20
=
0.4503 MPa at the temperature of 150°C.Figure 4
Stress strain curve.


## 2.1. Geometry and Meshing

The design parameters of the preform are shown in Figure1(a). The geometry is controlled by the three independent parameters: e, α, and H. Assume that the deformation of the preform occuring under constant volume (V
=
V
0) reduces the number of parameters from 3 to 2 as will be shown below. The FE analysis is performed using six-noded axisymmetrical elements with a full integration scheme. The preform is designed according to the pairs (
α
,
e
), Table 1. The mesh size is chosen in such a way that subsequent refinement does not affect the stress-strain behavior much. A typical FE mesh is shown in Figure 1(b).Table 1
Number of elements in the FE model.
α
/
e
3.6
3.8
4
4.2
2.1
4221/1872
4005/1776
3825/1696
3609/1600
2.3
—
4401/1952
4149/1840
3879/1720
2.5
—
—
—
4329/1920Geometry of the preform.
(a)
Geometrical parameterization
(b)
Typical mesh (zoom)The volume of the preform is calculated as the sum of a semisphere and a conical cylinder as shown in Figure1. After some algebraic manipulations, the expression of the volume of the preform is given by the three terms A, B, and S:
(1)
V
=
A
+
S
-
B
with
(2)
A
=
π
3
H
′
(
R
1
2
+
R
2
+
R
1
R
)
,
B
=
π
3
H
′
(
R
2
2
+
r
2
+
R
2
r
)
,
S
=
4
π
6
(
R
3
-
r
3
)
.
A, B, and S depend on the geometrical characteristics of the preform:
(3)
R
=
D
-
(
H
′
-
K
)
tan
α
,
r
=
(
D
-
e
)
-
(
H
′
-
K
)
tan
2
,
R
1
=
(
D
+
K
tan
α
)
,
R
2
=
(
D
-
e
+
K
tan
2
)
,
H
=
R
+
H
′
.
It is assumed that the deformation of the preform occurs under constant volume during the mold flow process (V
=
V
0). This leads to a third degree equation in e:
(4)
V
=
V
0
⟹
e
3
+
a
e
2
+
b
e
+
c
=
0
,
where
(5)
a
=
-
3
(
D
-
(
H
′
-
K
)
tan
2
)
-
3
2
H
′
,
b
=
3
(
D
-
(
H
′
-
K
)
tan
2
)
2
+
3
2
H
′
(
2
D
+
(
2
K
-
H
′
)
tan
2
)
,
c
=
3
2
π
(
A
-
V
0
)
+
1
π
(
R
3
-
(
D
-
(
H
′
-
K
)
tan
2
)
3
)
-
1
2
(
(
D
+
K
tan
2
)
2
+
(
D
-
(
H
′
-
K
)
tan
2
)
2
1
2
+
(
D
+
K
tan
2
)
(
D
-
(
H
′
-
K
)
tan
2
)
1
2
)
.
Analytical solution can easily be obtained for the third degree (4). Two cases are distinguished depending on the value of the determinant Δ
=
q
2
+
4
p
3
/
27 (q
=
2
a
3
/
27
-
a
b
/
3
+
c and p
=
b
-
a
2
/
3). For the studied problem Δ is negative, so that 3 real zeros verify (4). Each solution k (k
=
0, k
=
1, and k
=
3) is given by the means of polar coordinates:
(6)
e
=
2
r
1
/
3
cos
⁡
(
θ
3
+
2
k
π
3
)
with r
=
0.5
q
2
-
Δ and θ
=
arctan
(
-
-
Δ
/
q
). Only one solution lies on the allowed range. The diagram in Figure 2 can be divided into regions. The divisions between these regions are bounded by the curve with the equation H
=
H
max
⁡. The feasible region is above this curve, whereas the infeasible region below the curve leads to preform geometry in which the height H is greater than the height of the mold.Figure 2
Limit curve.


## 2.2. Boundary Conditions and Loading

The following boundary conditions are considered for the FE analysis.(i)
The “container” is fully constrained and the preform could move only along the axial direction.
(ii)
A displacement along thez-axis (Figure 3) to the bottom of the mold is applied first to the preform. This displacement is preconized by the manufacturer. Accordingly, the preform moves down to be in contact with the mold.
(iii)
A pressure of 0.3 MPa was applied in the second step to the inside of the preform with a gradual increase.
(iv)
A friction law, corresponding to Coulomb law with a constant friction coefficientμ
=
0.1, is considered between the preform and the “container.” It is worth noting that some simulations (see Section 3) in which μ was taken in the range of 0.1 to 0.5 have shown that the results were not affected by the value of the friction coefficient.Figure 3
Loading.


## 2.3. Material Behavior Law

To describe the mechanical behavior of the PMMA, hyperelastic and viscoelastic models are the two commonly used approaches. The emphasis of current study is a phenomenological hyperelastic approach due to its better suitability for the PMMA which is isotropic, highly elastic, and incompressible for the forming conditions [17]. It is already established that the Mooney-Rivlin model is accurate in simulating polymer sheet forming; see, for instance, [17, 18]. It should be noted that this approach fails to represent the speed influence. The material behavior is described with a specific free energy density W which is a function of strain and temperature. For isotropic behavior, W depends on the invariants of the transformation tensor F which is given by
(7)
d
x
=
F
d
X
,
B
=
F
T
F
,
where d
X is the elementary vector in the initial configuration, d
x is the transformation of d
X in the actual configuration, F is the transformation gradient tensor, and B is the Cauchy-Green Gauche tensor. The stresses can be obtained as the partial derivatives of the strain energy density function with respect to invariants I
1, I
2, and I
3. For isotropic materials, W is a linear combination of the invariants I
1, I
2, and I
3 of the strain tensor W
=
W
(
I
1
,
I
2
,
I
3
), where I
1
=
trace
B, I
2
=
(
1
/
2
)
(
(
trace
B
2
)
+
(
trace
B
)
2
), and I
3
=
det
⁡
(
B
). These invariants can be expressed with the elongations (λ
i, i
=
1,2, and 3) as follows:
(8)
I
1
=
λ
1
2
+
λ
2
2
+
λ
3
2
,
I
2
=
λ
1
2
λ
2
2
+
λ
1
2
λ
3
2
+
λ
2
2
λ
3
2
,
I
3
=
λ
1
2
λ
2
2
λ
3
2
.
The principal stretch ratio λ is the ratio of the gauge lengths after and before the test (λ
=
L
/
L
0). The hyperelastic theory (Mooney-Rivlin model) is adopted in this study for the simulation of the nonlinear material behavior. The strain energy density W is expressed by Mooney-Rivlin model as follows:
(9)
W
=
C
10
(
I
1
-
3
)
+
C
20
(
I
1
-
3
)
,
where C
10 and C
20 are two material parameters. For the case of incompressible material under uniaxial elongation (λ
1
=
λ, λ
2
=
λ
3
=
1
/
λ), the true Cauchy stress differences can be calculated as
(10)
σ
11
-
σ
33
=
2
C
10
(
λ
2
-
1
λ
)
+
2
C
20
(
λ
-
1
λ
2
)
,
σ
22
-
σ
33
=
0
.
In the case of simple tension, σ
22
=
σ
33
=
0, the uniaxial behavior law writes
(11)
σ
11
=
(
2
C
10
+
2
C
20
λ
)
(
λ
2
-
1
λ
)
.
The Mooney-Rivlin parameters of PMMA material are taken from the literature. One element test subjected to monotonic tensile loading is carried out in order to check the accuracy of the material parameters in simulating the experimental stress/elongation curve (Figure 4). It can be seen that the experimental data are well described by the Mooney-Rivlin model for the parameters values of C
10
=
-
0.0303 MPa and C
20
=
0.4503 MPa at the temperature of 150°C.Figure 4
Stress strain curve.


## 3. Results and Discussions

In this section, the FE results are discussed in terms of thickness distribution. Since the friction behaviour needed for the simulation of the stretch blow molding process is unknown, parametric study related to the effect of friction coefficient between the mold and the preform is preliminarily conducted. The prediction using three different friction coefficients, namely, 0.0 (frictionless), 0.1, and 0.5, is investigated. Figure5 shows that the variation of the friction coefficient does not affect the bubble thickness. Further simulations will be carried out with friction coefficients of 0.1. Figure 6 shows the distribution of the final product thickness versus the ordinate-abscissa for several geometries of the preform. Each geometry is characterized by a given pair (
α
,
e
) (given in Table 1). The thicknesses are in agreement with the results obtained by Cosson et al. [7] and Yang et al. [11]. The simulation results are superimposed to the measured data of an existing bubble. Figure 6 shows how the scattering parameters (e and α) can greatly improve the thickness uniformity. The standard deviation of the thickness is reduced to about 95% compared to the existing bubble. For a better evaluation of the effect of the geometrical parameters e and α of the preform on the final results, a contour plot of the standard deviation of the thickness is shown in Figure 7. Classically, an appropriate cost function is mathematically formulated in which the objective function depends on the standard deviation of the computed thicknesses. In our work, the number of driving parameters is reduced to 2 (i.e., e and α). It is then possible, thanks to the contour plot of the standard deviation, to optimize manually the thickness distribution. It can be seen that the thickness is widely spread for high values of e and low values of α. Low standard deviation showing that the thickness is clustered closely around a mean value is noticed for high values of e and low values of α. The increase of the preform thickness improves the thickness distribution of the bubble and decreases the height. Accordingly, the stretching increases and improves the mechanical properties of material of the finished product as demonstrated by Cosson et al. [8]. The “optimal” distribution of the thickness should be obtained for α
≃
2.1°C and e
≃
4.2 mm. It gives a light bubble with a uniform distribution of the thickness in addition to a better mechanical quality [7, 8].Figure 5
Effect of the friction coefficient.Figure 6
Bubble thickness evolution ordinate-abscissa.Figure 7
Contour plot of the standard deviation of the thickness (α versus e).


## 4. Conclusion

In this study, a finite element simulation is carried out to simulate the stretch blow molding. The geometrical parameters of the preform, leading to a homogeneous thickness distribution of the electric bubbles, are predicted. The simulation results show an improvement compared to an existing bubble in terms of thickness distribution and mechanical properties. Several possible extensions can be performed in the near future:(i)
the determination of the blow injection mold geometry allowing the realization of the electric bubbles,
(ii)
the choice of the mold cavity disposition, the alimentation system, and the cooling mechanism,
(iii)
the simulation of the process using the designed mold to optimise the implementation parameters such as cooling time, filling time, and closing force.



---

*Source: 175743-2014-08-27.xml*