# Simulation of Free Blowing of Polyethylene terephthalate Using a Thermodynamic Induced Crystallisation model

Luc Chevalier, Yann Marco, Guillaume Racineux, Samia Hmida Maamar

► To cite this version:

Luc Chevalier, Yann Marco, Guillaume Racineux, Samia Hmida Maamar. Simulation of Free Blowing of Polyethylene terephthalate Using a Thermodynamic Induced Crystallisation model. International Journal of Forming Processes, 2006, 9 (1), pp.29-59. 10.3166/ijfp.9.29-59. hal-05252765

HAL Id: hal-05252765

https://hal.science/hal-05252765v1

Submitted on 13 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Simulation of Free Blowing of Polyethylene terephthalate Using a Thermodynamic Induced Crystallisation model

## Free blowing simulation

L.Chevalier* - S.Maamar* - Y.Marco** - G.Racineux***

* LaM (EA 2545) - Marne la Vallée University, 5 boulevard Descartes, Champs sur Marne,F-77454 Marne la Vallée Cedex 2 France
luc.chevalier@univ-mlv.fr

** IUT de Saint-Brieuc 18 rue H. Wallon 22000 Saint-Brieuc France

** LEMA (CNRS FRE 2481) - Versailles Saint Quentin University, 45 avenue des Etats-Unis, 78035 Versailles France

RÉSUMÉ. L'évolution de la microstructure du poly(ethylene terephthalate) durant l'opération de soufflage a une influence importante sur les modifications de propriétés mécaniques (Chevalier, 1999). Pour mener à bien la simulation du procédé de soufflage par biororientation, il est nécessaire de prendre en compte ces évolutions. Dans une première partie, nous présentons le modèle de comportement (Poitou, 2001) dans le cas simple de viscosité en loi-puissance. Le modèle est identifié sur la base d'essais de traction uniaxiale et biaxiale régulés thermiquement. Dans la seconde partie, le problème axisymétrique de soufflage de membrane est résolu par éléments finis. Les simulations sont comparées aux préformes soufflées avec un succès encourageant compte tenu de la simplicité de la modélisation.

**ABSTRACT.** Microstructure evolution in Poly(Ethylene Terephthalate) during a blowing process has an important influence on final mechanical property (Chevalier, 1999). In order to simulate blow-moulding process it is necessary to take into account this evolution. In a first part we present a thermodynamic model in the case of a simple dependence of the viscosity under the crystallisation ratio (Poitou, 2001). The model is identified on uniaxial and biaxial tension tests results. In a second part, the axisymmetric problem of membrane inflation is solved using a classical finite element technique. The simulated results are compared with experimental results on free blown preforms and the comparison is quite accurate considering the elementary form of the induced crystallisation model.

MOTS-CLÉS: Soufflage, Simulation numérique, Cristallisation induite.

KEYWORDS: Blow Moulding, Numerical Simulation, Induced Crystallisation
# 1. Introduction

Stretching PolyEthylene Terephthalate (PET for short) at high strain rate over the glass transition temperature has a benefit effect on the material strength. In a previous paper (Chevalier et al., 1999) we presented results on the influence of stretch blow moulding parameters on the mechanical and microstructural properties of the final product, especially on the crystallinity ratio induced by stretching and the increase of Young modulus related to this crystallinity. The goal of the global study in which this work takes place, is the simulation of the PET morphology evolution during the blow moulding of bottles. Poitou and al. (Poitou, 2001) have proposed a thermodynamic model for induced crystallisation coupled to mechanical behaviour. This model has been validated for shear loading and for small strain rates tension tests (Poitou et al., 2003).

In the first part of this study, we show that the viscous form of this model fails to reproduce experimental results obtained at high strain rate in tension test of PET. We propose an identification procedure that shows that for a monotonous dependence of viscosity with respect to strain, the prediction of the model is not compatible with experimental results. This identification procedure is carried out considering uniaxial and biaxial results of tension tests managed by Marco (Marco et al., 2002). We propose a more simple form of the strain hardening effect and this form is used in a second part to simulate free blowing of a preform.

In this part, we present the basics of the kinematics description and the finite element type used to model the free blowing problem. Influence of the strain hardening effect is highlighted on the classical cylinder blowing problem which is also used to compare the finite element formulation with the analytical solution when viscosity is constant. Finally we show the effect of the strain induced crystallisation on shape evolution of the preform during free blowing.

# 2. Modelling Induced crystallisation in P.E.T.

The basics of the viscous form of the thermodynamic strain induced crystallisation model are briefly presented in this section. Uniaxial and biaxial tension tests results are presented and identification of the induced crystallisation model is discussed.

## 2.1. Basics of the strain induced crystallisation model

Thermodynamic study of material behaviour consists in the choice of accurate macroscopic variables to represent state of material. As presented in (Poitou et al.,
2001) we prefer, rather than the crystallinity ratio $\alpha$, the state variable $f$ called fictive crystallinity ratio and defined by Eq.1.

$$ \alpha = \frac{X_c}{X_c^\infty} = 1 - \exp(-f) \quad [1] $$

Dual variable to fictive crystallinity ratio $f$ is noted $\Lambda$. The expressions of the free energy and the pseudo-potential of dissipation give the behaviour and evolution laws by Eq.2.

$$ \Lambda = \frac{\partial \varphi}{\partial f} \quad \text{and} \quad \dot{f} = -\frac{\partial \varphi^*}{\partial \Lambda}; \dot{\varphi} = \frac{\partial \varphi^*}{\partial \Sigma} \quad [2] $$

Derivation of the free energy gives the relation between $f$ and its dual variable:

$$ \Lambda = \beta(T) \exp(-f) \quad \text{with:} \quad \beta(T) \approx LX_c^\infty \left( \frac{T}{T_m} - 1 \right) \quad [3] $$

$T$ is the current temperature and $L$ the latent heat of crystallisation. It worth noting that if we mix Eq.1 and Eq.3, the dual variable $\Lambda$ is directly related to the reduced crystallinity ratio $\alpha$ by Eq.4.

$$ \Lambda = \beta(T)(1-\alpha) \quad [4] $$

One can also see from Eq.3 that during isothermal transformation the kinetic of the state variable $\Lambda$ is related to the evolution of fictive crystallinity ratio $f$ by Eq.5.

$$ \dot{f} = -\frac{\dot{\Lambda}}{\Lambda} \quad [5] $$

Since infra red heating is quick and quiescent crystallisation slow for P.E.T., we can neglect the quiescent part in both free energy and pseudo potential. Nevertheless it is possible to introduce mechanical effects if viscosity $\eta$ depends on crystallinity ratio as in pseudo potential expression of Eq.6 for the linear case.

$$ \varphi^* = \frac{1}{4\eta(\Lambda)} \mathrm{Tr}(\Sigma^2) \quad [6] $$

It comes that second invariant of the extra-stress tensor is an engine for the crystallisation kinetic.

$$ \dot{f} = -\frac{\partial \varphi^*}{\partial \Lambda} = -\frac{\eta'(\Lambda)}{4\eta^2(\Lambda)} \mathrm{Tr}(\Sigma^2) \quad [7] $$

It is usually admitted that a constant viscosity does not give an accurate representation of the strain rate effect of polymer behaviour. We propose to consider a non linear form of the PET behaviour given by Eq.8 for the pseudo potential of dissipation. This expression gives a power law for viscosity.
$$ \varphi^* = \frac{m}{\frac{m+1}{2^{2m}} K^{\frac{m}{m+1}} (m+1)} \left[ \mathrm{Tr}(\Sigma^2) \right]^{\frac{m+1}{2m}} \quad [8] $$

By derivation we obtain the viscous part of the strain rate tensor which leads to:

$$ \Sigma = 2K \bar{\gamma}^{m-1} \dot{\bar{\epsilon}} \quad \text{with:} \quad \bar{\gamma} = \sqrt{2Tr(\dot{\bar{\epsilon}}^2)} \quad [9] $$

Slopes in the Log-log chart of Fig.3 are identical for all three curves so that exponent “m” can be chosen as a constant value. We assume that only “K” depends on crystallisation ratio. In that case, we obtain the fictive ratio evolution law by derivation of the pseudo-potential:

$$ \dot{j} = - \frac{\partial \varphi^*}{\partial \Lambda} = - \frac{K'}{\frac{m+1}{2^{2m}} K^{\frac{m}{m+1}} (m+1)} \left[ \mathrm{Tr}(\Sigma^2) \right]^{\frac{m+1}{2m}} \quad [10] $$

Where K' is the derivative of K with respect of associated variable Λ. Finally, the global problem consists to solve (i) the incompressibility condition; (ii) the strain rate-velocity relation; (iii) the equilibrium equation; (iv) the viscoelastic behaviour law and (v) the kinetic of crystallisation relation.

$$ \begin{gathered} \mathrm{div}(\underline{V}) = 0 \\ \dot{\underline{\epsilon}} = \frac{1}{2} \{ \underline{\mathrm{grad}} \underline{V} + \underline{\mathrm{grad}}^T \underline{V} \} \\ \mathrm{div} \underline{\Sigma} = - \underline{\mathrm{grad}} p \end{gathered} \quad [11] $$

$$ \begin{gathered} \underline{\Sigma} = 2K(\Lambda) \bar{\gamma}^{m-1} \dot{\underline{\epsilon}} \quad \text{with:} \quad \bar{\gamma} = \sqrt{2Tr(\dot{\underline{\epsilon}}^2)} \\ \frac{\dot{\Lambda}}{\Lambda} = \frac{K'(\Lambda)}{\frac{m+1}{2^{2m}} K(\Lambda)^{\frac{m}{m+1}} (m+1)} \left\{ \mathrm{Tr}(\underline{\Sigma}^2) \right\}^{\frac{m+1}{2m}} \end{gathered} $$

The identification problem to solve first is to determine exponent “m” and K(Λ) expression from experimental data. Relation between Λ and α are known since β(T) can be easily evaluated for PET. The identification problem is examined in the following section.

## 2.2. Tension tests on P.E.T.

In this section, we briefly present experimental results from uniaxial and biaxial elongation tests (Fig.1) carried out on injected specimen of PET during the Thesis work of Marco (Marco, 2003). The uniaxial and biaxial tension tests results are used to identify the strain induced crystallisation model.
<figure><img src="image_2.png" /></figure>

**Figure 1.** Infra red heating apparatus during uniaxial (a) and biaxial (b) tension tests. One can see infra red heating system on Fig.1a and the displacement field represented by vectors on Fig.1b. This displacement field is obtained by using a digital image analysis (Chevalier et al., 2001).

### 2.2.1 Experimental uniaxial tension tests

Specimens are injected using PET 9921W from Eastmann and following their recommendations. To minimise the effect of quiescent crystallisation, specimens are quickly heated with infra-red lamps before the test and temperature is regulated during the test. In order to analyse these effects on mechanical and microstructural properties we carried out tension tests at controlled temperature and velocity.

Tension tests have been carried out at different temperature (from 80°C up to 105°C) for a wide range of tension velocity (from 0.01 s⁻¹ to 1 s⁻¹). Each specimen is stretched from 50 mm up 200 mm (λ = 4). Typical results obtained for 90°C at the highest velocity available of 100 mm/s is presented on Fig.2a. Tension test is repeated and mean behaviour curve as well as lower and upper bound are plotted. One can see that dispersion occurs essentially during the strain hardening part (since 2.5 and 3 elongation level) of the test and this dispersion is important (about 20% of the mean stress value).

Figure 2b summarises all mean curves for different velocity at 90°C temperature. It worth noting that strain hardening always began between 2.5 and 3 elongation level. In the following we will focus on the influence of velocity on the behaviour of PET for this specific temperature. Considering the viscous aspect of the behaviour one can be interested to plot the apparent viscosity that can be calculated from Cauchy stress and strain rate by Eq.12.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Tension tests results on PET at T=90°C: (a) dispersion of PET behaviour for 100 mm/s tension speed, mean value of relative dispersion 19% (b) influence of tension speed on PET behaviour, strain hardening increases with tension speed.</figcaption></figure>

$$ \eta = \frac{\sigma}{3\dot{\epsilon}} \qquad [12] $$

During a tension test, strain rate decreases with respect to elongation $\lambda$ following Eq.13.

$$ \dot{\epsilon} = \frac{\dot{\epsilon}_0}{\lambda} \qquad [13] $$

$\dot{\epsilon}_0$ is the initial strain rate equal to $V/L_0$, V is the tension velocity and $L_0$ the initial length of the specimen. Data of Fig.2b can be used to estimate the viscosity $\eta$ as plotted on Fig.3. One can see that: (i) no master curve is highlighted so that
viscosity does not only depend on strain rate like it is usually admitted for injection simulation for example ; (ii) in this log-log chart a power law could be identified for each velocity but velocity is not a material parameter. This means that viscosity is also affected by another parameter (crystallinity or strain or macromolecular orientation for example). Considering that viscosity follows a power law, one can calculate the tension stress as a function of strain rate by Eq.14.

$$ \sigma(t) = 3^{\frac{m+1}{2}} K \dot{\epsilon}^m(t) \Rightarrow K = \frac{\sigma}{3^{\frac{m+1}{2}} \dot{\epsilon}^m} \quad [14] $$

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Apparent viscosity versus strain rate. No master curve appears, the classical decreasing of viscosity depends on more than strain rates</figcaption></figure>

Equation 14 allows function K to be plotted with respect to strain ε for example (the best would be to plot K versus Λ or α) for each value of exponent “m”. Figure 3 showed that no master curve arises for the viscosity “η” versus strain rate but for specific values of “m” the parameter K appears to be independent of tension velocity V and only function of logarithmic strain ε (ε = Lnλ) as shown on Fig.4b. The optimisation of the “m” parameter value is managed by summing the difference between all experimental curves and a minimal dispersion is obtained for “m” equal to 0.6 as shown on Fig.4a. Thus we can identify a simple form for the K(ε) function following the relation given by Eq.15 with parameters:
$A = 0.916 \ 10^8$; $B = -2.986 \ 10^8$; $C = 3.705 \ 10^8$; $D = -2.197 \ 10^8$; $E = 0.642 \ 10^8$; $F = -0.044 \ 10^8$; $G = 0.0014 \ 10^8$.

$$K(\varepsilon) = A\varepsilon^6 + B\varepsilon^5 + C\varepsilon^4 + D\varepsilon^3 + E\varepsilon^2 + F\varepsilon + G \quad [15]$$

<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Minimization of differences between $K(\varepsilon)$ functions. An optimal value is obtained for $m=0.6$ (b) $K$ functions represented (circles) the difference between the three functions is lower than the dispersion of each. The bold curve is a polynomial identification of the mean value of $K$ function.</figcaption></figure>

Measured crystallinity ratios after the tension test for each specimen are given in table 1. It appears that tension velocity has a low effect on final crystallisation for this temperature.

<table><caption>**Table 1.** Induced crystallisation in PET specimen after tension test at 90°C</caption><tr><td>Tension velocity V (mm/s)</td><td>15</td><td>50</td><td>100</td></tr><tr><td>Crystallinity ratio Xc en (%)</td><td>35</td><td>37</td><td>39</td></tr><tr><td>α</td><td>0,7</td><td>0,74</td><td>0,78</td></tr></table>

### 2.2.2. Experimental biaxial tension tests

Biaxial tension tests have been managed by Marco (Marco, 2003) on biaxial testing machine Astree. In the equi-biaxial tension test, the square specimen has been simultaneously stretched in both directions X and Y. The principal elongation $\lambda_x$ and $\lambda_y$ in the central region of the specimen are equal and assure the equi-biaxial
state in this zone. One can also observe that the strain is quasi constant along AB (Fig.5a). As shown in (Chevalier et al., 2002), if the stress components are quasi uniform in this section AB we can estimate the biaxial Cauchy stress $\sigma_{BT}$ from the loading force F:

$$ \sigma_{BT} = \sigma_{nn} = \frac{\sqrt{2F}}{eL} \quad [16] $$

where “e” is the current specimen thickness and “L” the current length between A and B. Strain hardening effect occurs at the same elongation level than for uniaxial tests.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Biaxial tension tests results, biaxial stress is plotted version elongation in direction X or Y. Strain hardening effect occurs at same elongation level than for uniaxial tests.</figcaption></figure>

Using the classical definition of equivalent strain will make identification of a unique K(ε) curve impossible. In that case, the equivalent strain would be equal to “2Lnλ” which would put the strain hardening beginning in elongation level from 4 to 5 instead of 2.5 to 3. We prefer to consider as an equivalent strain definition the expression given by Eq.17.

$$ \underline{\epsilon} = \text{Sup}_i \epsilon_i \text{ for } i = 1, 2, 3 \quad [17] $$

Where $\epsilon_i$ are the three principal strains. From this biaxial stress $\sigma_{BT}$, one can obtain the K(ε) function by Eq.18.

$$ \sigma_{BT}(t) = 2^m 3^{\frac{m+1}{2}} K \dot{\epsilon}^m(t) \Rightarrow K = \frac{\sigma_{BT}}{2^m 3^{\frac{m+1}{2}} \dot{\epsilon}^m} \quad [18] $$
For a better identification of K(ε) function, the same procedure previously described has been carried out using the six K(ε) functions -3 from uniaxial tension test and 3 more from biaxial tension tests-. One can see that the optimum “m” value is still equal to 0.6 (Fig.6a) and that dispersion on K function is less than 1% for this value of exponent “m”. One can also use this optimal value to fit a polynomial form for K(ε). All parameters A, B, C etc... are given below and are very close to the values previously obtained from purely uniaxial tension tests.

$$K(\varepsilon) = A\varepsilon^6 + B\varepsilon^5 + C\varepsilon^4 + D\varepsilon^3 + E\varepsilon^2 + F\varepsilon + G \quad [19]$$

$$\text{with: } A = 0.854 \ 10^8; B = -2.745 \ 10^8; C = 3.399 \ 10^8; D = -2.012 \ 10^8;$$

$$E = 0.581 \ 10^8; F = -0.0442 \ 10^8; G = 0.0013 \ 10^8$$

<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Minimization of differences between K(ε) functions in the biaxial tension case. An optimal value is obtained for m=0.6 which is closer of the optimal value of the uniaxial case rather than the biaxial one. The global dispersion between all K function is equal to 0.83% (b) K functions for both uniaxial and biaxial cases are represented (circles). Differences between the six functions are small. The bold curve is a polynomial identification of the mean value of K function. (c) same than b with exponential function.</figcaption></figure>
Using a more specific form of K(ε) function, one can reduce the number of parameters needed. Figure 6b shows the identification done with an exponential form given by Eq.20.

$$K(\varepsilon) = K_0 \exp(A\varepsilon^B) \quad [20]$$

Considering the accurate prediction of biaxial behaviour from uniaxial data, we can conclude that the highest principal strain governs the strain hardening effect presented during elongation of PET. Sequential biaxial tests should be carried out to confirm this conclusion. Crystallinity ratio measured after the biaxial tension test are given in table 2. Once again it appears that velocity has not a great influence on crystallisation, but one can also observe that Xc values reached are quasi identical to uniaxial tests.

<table><caption>**Table 2.** Induced crystallisation in PET specimen after biaxial tension test at the temperature of 90°C</caption><tr><td>Tension velocity V (mm/s)</td><td>15</td><td>30</td><td>100</td></tr><tr><td>Crystallinity ratio Xc en (%)</td><td>36</td><td>37</td><td>42</td></tr><tr><td>α</td><td>0,72</td><td>0,74</td><td>0,84</td></tr></table>

### 2.2.3. Discussion on model identification

In this section, we will test the induced crystallinity model presented above and we will discuss the difficulties to prove that the crystallisation kinetic is responsible for the strain hardening effect.

First let's examine the crystallisation evolution law of the strain induced crystallisation model. Replacing the associated variable Λ by its expression with reduced crystallinity ratio α yields to Eq.21 to govern crystallisation kinetic.

$$\dot{\alpha} = \frac{3^{\frac{1+m}{2}} (1-\alpha) \varepsilon^{1+m}}{(m+1)\beta(T)} \frac{\partial K}{\partial \alpha} \quad [21]$$

The first remark is that crystallisation rate depends on (1-α). When α reaches 100% the parameter K becomes constant. Because of the decreasing form of the strain rate in a tension test, stress will also decrease as shown on one of the curves of Fig.7a. This evolution is not observed during tests so this means that crystallisation must be far from saturation. One can also notice that the strain hardening effect does not occur at the same elongation level for each tension velocity.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Typical behaviour prediction using the thermodynamic strain induced crystallisation model (b) Probabilistic density function of initial reduced crystallinity ratio (c) Dispersion on behaviour curves due to initial microstructure of specimen.</figcaption></figure>

Considering an analogous form for K(α) than the one identified for K(ε) one can try to fit the best triplet (K₀, Aα, B) for the expression of Eq.22.

$$K(\alpha) = K_0 \exp(A_{\alpha} \alpha^B) \Rightarrow \frac{\partial K}{\partial \alpha} = K_0 A_{\alpha} B \alpha^{B-1} \exp(A_{\alpha} \alpha^B) \quad [22]$$

Since crystallisation rate depends on ∂K/∂α, one can show that crystallisation rate is null if initial crystallinity ratio is null. In that case, the prediction of the strain induced model is identical to classical viscous model and no strain hardening effect arises. In fact, initial crystallinity ratio is not strictly equal to zero in injected specimen or preform. Initial value of Xc is about 4% and dispersion occurs from a specimen to another. It worth noting that this dispersion has an important effect in the model prediction. For example, if the reduced initial crystallinity follows a probabilistic law described by the probability density function of Fig.7b, one can show that the model will describe behaviour by the curves of Fig.7c. This dispersion is too important in regards of the experimental data. The model is very sensitive to
initial values of parameters and leads to numerical instabilities when used to simulate forming process.

In situ measurement managed at low strain rate (Chaari et al. 2002) shows that $\alpha$ can be related to equivalent strain $\varepsilon$ by a linear function. Such a form has already been proposed in (Vigny et al. 1999) and one can determine the slope by examining final values of $X_c$ measured after the tension test.

$$ \alpha = k \varepsilon \text{ with } k = 0.5 \qquad [23] $$

From expression of Eq.23, one can use Eq.20 to express $K(\alpha)$ in the form of Eq.22 with: $A_\alpha = A/k^B$. The derivative form of Eq.22 can be compared to the expression issued of Eq.21.

$$ \frac{\partial K}{\partial \alpha} = \frac{(m+1)\beta(T)}{3^{\frac{1+m}{2}} (1-\alpha) \dot{\varepsilon}^{1+m}} \dot{\alpha} \Rightarrow \frac{\partial K}{\partial \alpha} = \frac{k(m+1)\beta(T)}{3^{\frac{1+m}{2}} (1-\alpha) \dot{\varepsilon}_0^m} \exp\left(\frac{m\alpha}{k}\right) \qquad [24] $$

One can see that this expression depend on initial strain rate (i.e, $V/L_0$). If we can manage to fit the two expressions for one velocity, it won't fit for other velocities. This confirms the observation made on Fig.7a.

### 2.2.2.4. Weak form of the strain induced crystallisation model

A weaker form of the model is proposed to manage accurate simulation. The crystallisation kinetic is only governed by the maximal principal logarithmic strain. This assumption is only legal for a range of velocity, tests done at very low strain rate lead to lower crystallinity ratio (Marco, 2003). At very high strain rate, the strain hardening effect only depends on orientation of macromolecular chains as shown in (Mahendrasingam et al. 1999). Eq.21 is replaced by a much simple form:

$$ \dot{\alpha} = k \dot{\varepsilon} \Rightarrow \alpha = k Ln(\sup_i(\lambda_i)) \qquad [25] $$

$\lambda_i$ is the highest principal elongation. The flow problem is governed by the same relations than in Eq.11, where $K(\Lambda)$ is replaced by $K(\alpha)$ of Eq.23. Applied to uniaxial and biaxial tension test, this weak form of the Poitou's model gives the curves of Fig.8. It is an accurate representation of behaviour curves of Fig.2 and Fig.5.

In the following section, we will apply this model to the simulation of: (i) tube stretching and blowing; (ii) preform blowing. A membrane theory will first be developed and compared to a classical 3D modelling. Since geometrical evolution is similar for both modelling we will implement this theory in a finite element code for the blow moulding process simulation.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Simulation of tension test with the strain induced crystallisation model (a) Uniaxial tension test (b) Biaxial tension test.</figcaption></figure>

### 3. Free blowing of P.E.T. preforms

#### 3.1. Theoretical topics

Considering the low thickness ratio ($e/r$) of the blown bottle, we neglect bending effects and build a thin body tension theory. The problem (Fig.9) is axisymmetric so that we work in the meridian plane $(O, \vec{e}_r, \vec{e}_z)$. In this plane, a point $m$ on the middle surface of the bottle is described, at time $t$, by the coordinates $(r(s), \theta, z(s))$ where $s$ is the curve abscissa of $m$ on the meridian line.

$$ \overline{Om}(s, \theta) = r(s)\vec{e}_r + z(s)\vec{e}_z \quad [26] $$

In the local base $(\vec{n}, \vec{t}, \vec{s})$ the velocity at point $m$ writes:

$$ \vec{V}(m) = u(s)\vec{n} + w(s)\vec{s} \quad [27] $$

We suppose that tangential velocity (following $\vec{s}$) and normal velocity (following $\vec{n}$) are respectively uniform and linear on the bottle thickness. So, the velocity of a point $M$ at a distance $\xi$ of point $m$ (Fig.9) writes:

$$ \vec{V}(M) = (u(s) + \lambda(s)\xi)\vec{n} + w(s)\vec{s} = U\vec{e}_r + W\vec{e}_z + \lambda\xi\vec{n} \quad [28] $$

$\lambda$ function is obtained by the incompressibility condition.

Supposing that $\xi$ is small compare to radius $r$ and curvature radius $R_s$, the strain rate tensor writes:
$$ \dot{\underline{\underline{\epsilon}}} = \begin{pmatrix} \lambda & 0 & \frac{1}{2} \left( u' + \lambda' \xi - \frac{w}{R_s} \right) \\ 0 & \frac{(u + \lambda \xi) \cos \alpha - w \sin \alpha}{r} & 0 \\ \frac{1}{2} \left( u' + \lambda' \xi - \frac{w}{R_s} \right) & \frac{w'}{R_s} & \frac{(u + \lambda \xi)}{R_s} \end{pmatrix}_{(\dot{e}_n, \dot{e}_r, \dot{e}_s)} \quad [29] $$

with

$$ R_s = \frac{1}{z''r' - r''z'} \quad [30] $$

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Axisymmetric thin body tension theory. (a) mean surface description, (b) thickness description</figcaption></figure>

Primes and double primes indicate first and second derivatives with regard to curve abscissa. According to cinematic hypothesis, the incompressibility condition must be satisfied on thickness *e* which leads to:
$$ \langle \mathrm{div}(\vec{V}(M)) \rangle_e = \frac{1}{e} \int_{-e/2}^{e/2} \mathrm{div}(\vec{V}(M)) d\xi = \lambda + w' + \frac{u}{R_s} + \frac{u \cos \alpha - w \sin \alpha}{r} = 0 \quad [31] $$

### 3.2. Governing equations and variational formulation

The solution fields researched $(\vec{V}, p)$ are respectively element of the spaces $V_v$ and $V_p$ with:

$$ V_v = \{ \vec{V} \in H^1(\Omega) / \vec{V} = \vec{V}_d \text{ on } \partial\Omega_v \} \quad [32] $$

$$ V_p = L^2(\Omega) \quad [33] $$

We consider variations of solution field such as

$$ \vec{V}^* \in V_v^0 = \{ \vec{V} \in H^1(\Omega) / \vec{V} = \vec{0} \text{ on } \partial\Omega_v \} \quad [34] $$

The problem consist to find the velocity field $\vec{V}$, the thickness $e$ and the pressure field $p$ verifying:

- balance equations

$$ \int_{\Omega} \langle \mathrm{div}(\vec{V}) \rangle_e p^* dv = 0, \quad \forall p^* \in V_p \quad [35] $$

$$ \int_{\Omega} \mathrm{tr}[\underline{\underline{\Sigma \dot{\epsilon}^*}}] - \int_{\Omega} p \langle \mathrm{div}(\vec{V}^*) \rangle_e dv - \int_{\partial \Omega_p} \vec{F}^d \cdot \vec{V}^* dv = 0, \quad \forall \vec{V}^* \in V_v^0 \quad [36] $$

- constitutive equation

$$ \underline{\underline{\Sigma}} = 2\eta\dot{\underline{\underline{\epsilon}}}} \quad \text{with} \quad \eta = 2K(\alpha)\left(\sqrt{2\mathrm{tr}(\dot{\underline{\underline{\epsilon}}^2})}\right)^{m-1} \quad [37] $$

where

$$ K = K_0 \exp(A\bar{\epsilon}^B), \quad \bar{\epsilon} = \mathrm{Sup}(\mathrm{Ln}\lambda_i) \quad [38] $$

- boundary and initial conditions

$$ \vec{V} = \vec{V}_d, \quad \text{on } \partial\Omega_v \quad [39] $$

$$ \vec{F}_d = -\Delta p \vec{n} \quad \text{on } \partial\Omega_p \quad [40] $$

$$ e(s) = e_0(s) \quad \text{at } t=0 \quad [41] $$

Considering constitutive law (Eq. 37) and boundary conditions (Eq. 39-40), the equivalent variational formulation of balance equations writes:
Finding $(\vec{V}, p) \in (V_v \times V_p)$ such as $\forall (\vec{V}^*, p^*) \in (V_v^0 \times V_p)$:

$$ \int_{\Omega} \mathrm{tr} \left[ 2\eta \underline{\underline{\dot{\epsilon} \dot{\epsilon}^*}} \right] dv - \int_{\Omega} p \left\langle \mathrm{div} \left( \vec{V} \right) \right\rangle_e dv - \int_{\partial \Omega_p} \vec{F}^d \cdot \vec{V}^* dv - \int_{\Omega} \left\langle \mathrm{div} \left( \vec{V} \right) \right\rangle_e p^* dv = 0 \quad [42] $$

### 3.3. Tube inflation simulation

The analytical solution of the classical problem of cylinder inflation under membrane hypothesis and for a Newtonian behaviour is presented in the appendix with local equilibrium equations coming out of this theory. In this section, we examine the influence of the strain hardening of PET for same problem, using a power law viscosity with K function of $\epsilon$. Problem to solve writes:

$$ \left\{ \begin{aligned} \dot{R} &= \frac{\Delta P}{4\eta} \left( \frac{R}{e} + \frac{1}{2} \right) R - \frac{\dot{L}R}{2L} \\ \dot{e} &= -\frac{\Delta P}{4\eta} \left( \frac{R}{e} + \frac{1}{2} \right) e - \frac{\dot{L}e}{2L} \end{aligned} \right. \qquad [43] $$

with $\eta = K \dot{\gamma}^{m-1}$ where $\dot{\gamma} = \sqrt{2\mathrm{trace}(\underline{D}^2)}$

and $K = K_0 \exp(A \bar{\epsilon}^B)$ where $\bar{\epsilon} = \mathrm{Sup}(\mathrm{Ln} \lambda_i)$.

Because the problem is no longer linear, we use a Newton method to solve the flow problem at each time step:

$$ \left\{ \begin{aligned} \left( \begin{matrix} \dot{R} \\ \dot{e} \end{matrix} \right)_{t(n+1)} &= \left( \begin{matrix} \dot{R} \\ \dot{e} \end{matrix} \right)_{t(n)} - \left( \overline{\mathrm{grad}}(F_t(n)) \right)^{-1} F_t(n) \\ \dot{L}_{t(n+1)} &= \dot{L}_{t(n)} = W \end{aligned} \right. \qquad [44] $$

$$ \text{where} \quad \left\{ \begin{aligned} F_1(\dot{e}, \dot{R}) &= \dot{R} - \frac{\Delta P}{4\eta} \left( \frac{R}{e} + \frac{1}{2} \right) R + \frac{\dot{L}R}{2L} \\ F_2(\dot{e}, \dot{R}) &= \dot{e} + \frac{\Delta P}{4\eta} \left( \frac{R}{e} + \frac{1}{2} \right) e + \frac{\dot{L}e}{2L} \end{aligned} \right. \qquad [45] $$

and

$$ \overline{\mathrm{grad}}(F) = \begin{pmatrix} 1 + \frac{2(m-1)}{K} \frac{\Delta P}{4} \left( \frac{R}{e} + \frac{1}{2} \right) \frac{\dot{R}}{R} \dot{\gamma}^{-(m+1)} & \frac{2(m-1)}{K} \frac{\Delta P}{4} \left( \frac{R}{e} + \frac{1}{2} \right) \frac{R \dot{e}}{e^2} \dot{\gamma}^{-(m+1)} \\ -\frac{2(m-1)}{K} \frac{\Delta P}{4} \left( \frac{R}{e} + \frac{1}{2} \right) \frac{\dot{R} e}{R^2} \dot{\gamma}^{-(m+1)} & 1 - \frac{2(m-1)}{K} \frac{\Delta P}{4} \left( \frac{R}{e} + \frac{1}{2} \right) \frac{\dot{e}}{e} \dot{\gamma}^{-(m+1)} \end{pmatrix} \qquad [46] $$
We consider a cylinder with initial geometry given by: $R_0 = 10mm$, $L_0 = 45mm$, $e_0 = 4mm$ that can be representative of a PET preform. On Fig.10, the aspect ratio evolution versus time is plotted for both cases: classical power law viscosity ($K = K_0$) and power law viscosity coupled with induced crystallisation. In the first case, $f(t)$ decreases to zero and blowing presents an ultimate time as for Newtonian viscosity. This is no longer the case if induced crystallisation is taken into account. Radius continuously increases and seems to reach an asymptotic value.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** (a) Aspect ratio evolution during tube inflation with and without considering the stress hardening effect (b) radius and thickness evolution.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Aspect ratio evolution during tube inflation with stress hardening effect for different longitudinal speeds $W$ (b) radius and thickness evolution.</figcaption></figure>

During stretch-blow moulding process, the bottle (cylinder) is stretched along the longitudinal direction while blowing occurs. This case is represented on Fig.11 for different value of the elongation velocity $W$ (0, 1 and 10 m/s). The two first
cases are quite similar but when longitudinal elongation becomes high, the radius evolution is very slow because of the elongation and thickness quickly decreases. Then blowing is easier but since strain hardening increases and longitudinal elongation continues, then radius decreases.

These cases highlights the influence of the strain hardening effect on geometry evolution but realistic boundary condition can only be applied if we described accurately the preform geometry. This can be done using finite element simulation and is presented in the next section.

### 3.4. Finite element simulation

Since the problem is unsteady, we adopted an incremental approach in three steps:

In the first step, we assume that the thickness $e$ is known and we solve the flow problem $(\bar{V}, p)$. We use the finite element method with the velocity-pressure formulation (Eq. 42) and a Newton-Raphson algorithm to solve the non-linearity.

Then, the thickness is actualised using finite difference.

At the end of the time increment the mesh is convected.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Axisymmetric linear element</figcaption></figure>

The problem is axisymmetric and we use membrane elements so, the 3D problem can be reduced to a 1D axisymmetric problem. For geometric interpolation we use a linear element (Fig.12), then $\alpha = \text{cst}$ on an element and $R_s \to +\infty$. To avoid locking problem and to keep an homogenous interpolation degree between stress tensor terms, we need to use quadratic interpolation (SEG3) for the velocity and linear (SEG2) one for the pressure. The following numerical simulation are carried out from a simplified preform defined by a half spherical part and a cylindrical one. The mean radius of both the cylinder and the sphere is $r_0 = 10$ mm.
The cylinder length is $L_0 = 45$ mm. The global length ($L_0 + \pi r_0 / 2$) is divided in 200 finite elements.

First let's examine the case described in appendix, where Newtonian behaviour is chosen for the simulation. One can note (Fig.20) that the shape evolution of the blown preform yields to a spherical bottle and is in contradiction with the photograph of Fig.16. The asymptotic time evolution is the consequence of the lack of strength due to the thickness decreasing combined with the increasing of internal surface for the pressure. The evolution of the shape is due to the equilibrium of the preform, it can be shown that internal pressure generates an axial load which is not sufficient to elongate the preform. Consequently, the radius increases in the central region of the cylinder and a spherical shape arises. When this evolution has started, nothing can lead to bottle shape, only a strain hardening effect could stop geometrical evolution when strain in the circumferential direction is higher than in the longitudinal one.

Second, using parameters identified on tension tests where viscosity depends of crystallisation ratio, three cases are studied:

In the first case, initial thickness is uniform $e_0 = 4$ mm and internal pressure is constant $P = 10^6$ Pa. The simulation gives the evolution shown on Fig.13. One can note that the preform tend to stretch more than inflate, this is due to the effect of the strain induced crystallisation. The shape evolution of the blown preform is similar to the photograph of Fig.16. Radius evolution is the same that in the analytical case of the tube presented in the preceding paragraph.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** (a) Preform evolution, case 1: uniform initial thickness $e_0 = 4$ mm and constant internal pressure $P = 10^6$ Pa (b) Maximal radius and maximal length of the bottle versus time</figcaption></figure>
The radial blowing begins in the central region and the PET becomes more viscous so that deforming the bottom is easier than growing the radius. The shape evolution leads to a long bottle according to the experimental results shown on picture of Fig.16. The inflation of the shape becomes more and more difficult as viscosity increases and finally saturates.

The average crystallinity ratio grows up to 30% during simulation of this process. Density measured made on cut out samples from both ‘bottles’ leads to an average crystallinity ratio of 32% which is a quite accurate prediction. Considering the limitation due to the use of a membrane element and the hypothesis of constant temperature, it is not possible to be more quantitative in the comparison.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** (a) Preform evolution, case 2: uniform initial thickness $e_0 = 4$ mm and constant mole rate 10 moles/sec, $N_0 = 5.4 \ 10^{-3}$ moles (b) Maximal radius and maximal length of the bottle versus time (c) Internal pressure evolution versus time.</figcaption></figure>

Nevertheless, improvement in the geometrical modelling and loading can be done. Constant internal pressure is certainly not representative of the blowing process, as done in (Erchiqui et al., 2004), we used a perfect gas approximation for internal air behaviour. The law of perfect gas given by Eq.47, internal pressure is variable during free blowing and its evolution is shown on Fig.14c if initial thickness is uniform and mole ratio is constant.
$$PV = NRT$$

Where: “P” is internal pressure, “V” the volume occupied by the gas, “N” the mole number, “T” the temperature and “R” is the constant of perfect gas.

The simulation gives the evolution shown on Fig.14. One can note that the preform is stretched more quickly than in the first case, this is due to the effect of the variation internal pressure during free blowing.

In the final case (Fig.15), initial thickness is more representative of an industrial preform. The mole ratio is constant but initial thickness varies near the neck of the bottle and in the bottom zone. One can note that the evolution of the preform shape is different of the other cases: in spite of the internal pressure variation during free blowing, the preform is stretched less quickly than the preceding case and this is due to the effect of non uniform thickness.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** (a) Preform evolution, case 3: non uniform initial thickness $e_0$ = varying from 2 to 4 then to 3 mm and constant mole rate (idem Fig.14) (b) Maximal radius and maximal length of the bottle versus time (c) Internal pressure evolution versus time.</figcaption></figure>
Comparison between simulated blown shapes and real blown preforms shows that final shapes are well represented since strain hardening is taken into account. Further comparison between thickness distribution, crystallinity ratio distribution also shows good qualitative agreement. As it has already been mentioned, numerical simulation needs to be improved to make quantitative comparison possible.

<figure><img src="image_18.png" /><figcaption>**Figure 16.** *Preform and free blown samples*</figcaption></figure>

**4. Conclusions and perspectives**

Using experimental data from tension test made on P.E.T., identification of a weak form of a thermodynamic strain induced crystallisation model have been managed successfully. The behaviour chosen is simple and need few adjustable parameters ($K_0$ and $m$ for the power law and A and B for the dependence of viscosity wrt. crystallinity ratio).

A membrane theory has been presented and implemented in a finite element program in order to simulate free blowing of preform. This finite element code have been validated on the classical case of tube inflation. Results clearly show the influence of strain hardening (related to crystallinity) on the preform shape evolution.

As a perspective, a viscoelastic version of the model is implemented since this behaviour fits much better experimental data and the polymer behaviour. The kinetic law of crystallinity should be improved to take into account the strain rate effect on a wider range. To take into account the influence of temperature and radial elongation that may be different from inside to outside of the preform, numerical simulation will be managed using 2D axisymmetric elements. A free blowing apparatus is also developed in order to make more quantitative comparison between numerical simulation, blown preform shape and microstructure.
Acknowledgement

Sidel provided 9921W granular P.E.T. for tension specimen injection.
Specimens have been injected at LTVP ENSAM Paris.

References

Chaari F., Chaouche M., Chevalier L., Doucet J., "Real-time investigation of the strain-induced crystallization of poly(ethylene terephtalate) using synchrotron radiation : crystallization state evolution versus mechanical behavior", *6th European Rheology Conference* 2002, Erlangen (Allemagne)

Chevalier L., Calloch S., Hild F., Marco Y., "Digital correlation used to analyse the multi-axial behaviour of rubber-like materials", *Eur. J. Mech. A/solids*, 20, 169-187 (2001)

Chevalier L., Linhone C., Régnier G., "Induced Crystallinity During Stretch-blow Molding Process and its Influence on Mechanical Strength of Poly(Ethylene Terephthalate) Bottles", *Plastics Rubber & Composites, Processing & application*, 28, 8, p393-400 (1999)

Chevalier L., Marco Y., "Tools for Multiaxial Validation of Behaviour Laws Chosen for Modeling Hyper-elasticity of Rubber-like Materials", *Polymer Engineering & Science*, issue vol.42, n°2 (2002)

Erchiqui F., Gakwaya A., Diri D., Majid B., "Analysis of gas pressure during the thermoplastic membrane blowing using the dynamic finite element method", to be published in Int. J. Forming Proc. 2004

Mahendrasingam A., Martin C., Fuller W., Blundell D.J., MacKerron D. H., Oldman J., Harvie J. L., Riekel R. C., Engström P., "Effect of draw ratio and temperature on the strain-induced crystallisation of PET at fast draw rates", *Polymer* 40, 5553-5565 (1999)

Marco Y., « Caractérisation multi-axiale du comportement et de la microstructure d'un semi-cristallin : application au cas du PET », Thèse de doctorat de l'ENS de Cachan, Juin 2003

Marco Y., Chevalier L., Régnier G. and Poitou A., "Induced Crystallisation and orientation of Poly(ethylene-terephthalate) during Uniaxial and Biaxial elongation", *Macromolecular Symposia*, 185, p15-34, (2002).

Poitou A., Ammar A., "Polymer Crystallisation Induced by a Strain or Flow: a Thermodynamic approach", *C.R.Acad.Sci.* Paris, 329, série IIb, p5-11 (2001)

Poitou A., Ammar A., Marco Y., Chaouche M., Chevalier L., "Crystallisation induced by strain and flow of polymer : from molecular properties to macroscopic models", *Comp. Methods in Mechanics* 192, p3245-3264 (2003)

Vigny M., Aubert A., Hiver J-M., Aboulfaraj M., G'Sell C., "Constitutive viscoplastic behaviour of amorphous PET during plane-strain tensile stretching", *Polymer Eng. And Science*, 39 (12) 2366. (1999)
# Appendix: Local equations and tube inflation example

In order to compare the numerical approach, we performed the simulation of the inflation of a cylindrical tube. In case of Newtonian behaviour, the local equations can be solved analytically and the numerical solution can be compared with that reference. In this appendix, we first establish local relations of membrane theory and then study the tube inflation case.

## Stress definition and local relation

As shown previously, on the local base $(\vec{n}, \vec{s}, \vec{t})$, the velocity field writes:

$$ \vec{V}(M) = \{u(s) + \lambda(s)\xi\}\vec{n} + w(s)\vec{s} \quad [1] $$

Assuming $\xi$ is small in regard to radius $r$ or $R_s$, the strain rate tensor expression is:

$$ \dot{\mathbf{\epsilon}} = \begin{bmatrix} \lambda & 0 & \frac{1}{2}\left(u' + \lambda'\xi - \frac{w}{R_s}\right) \\ 0 & \frac{(u + \lambda\xi)\cos\alpha - w\sin\alpha}{r} & 0 \\ \frac{1}{2}\left(u' + \lambda'\xi - \frac{w}{R_s}\right) & 0 & w' + \frac{u + \lambda\xi}{R_s} \end{bmatrix}_{\vec{n}, \vec{t}, \vec{s}} \quad [2] $$

The incompressibility condition cannot be satisfied at all point M. It can be written over the thickness “e” which leads to:

$$ \langle \mathrm{div} \vec{V} \rangle = \lambda + \frac{u \cos \alpha - w \sin \alpha}{r} + w' + \frac{u}{R_s} = 0 \quad [3] $$

Using principle work principal, we can define the extra-stress tensor $\Sigma$:

$$ - \int_{\Omega} \mathrm{Tr}(\Sigma \dot{\mathbf{\epsilon}}^*) dV + \int_{\Omega} p \langle \mathrm{div} \vec{V}^* \rangle dV + \int_{\partial \Omega_f} \Delta P \vec{n} \cdot \vec{V}^* dA = 0 \quad [4] $$

A stress assumption is made to consider shear stress as neglecting (i.e., the extra-stress tensor is diagonal). Developing the different parts under the kinematics assumption, one can find the following expression:

$$ \mathrm{Tr}(\Sigma \dot{\mathbf{\epsilon}}^*) = \Sigma_n \lambda^* + \Sigma_t \frac{((u^* + \lambda \xi)\cos\alpha - w^*\sin\alpha)}{r} + \Sigma_s \left(\frac{dw^*}{ds} + \frac{u^* + \lambda \xi}{R_c}\right) \quad [5] $$

Integration over the thickness and over $\theta$ leads to the expression:

$$ \int_{\Omega} \mathrm{Tr}(\Sigma \dot{\mathbf{\epsilon}}^*) dV = 2\pi e \int_{s_1}^{s_2} \left\{ r \Sigma_n \lambda + \Sigma_t (u^* \cos\alpha - w^* \sin\alpha) + r \Sigma_s \left( \frac{dw^*}{ds} + \frac{u^*}{R_c} \right) \right\} ds \quad [6] $$

Integration by parts of the $dw^*/ds$ term finally leads to the following expression:
$$ \int_{\Omega} \mathrm{Tr}(\Sigma \cdot \dot{\epsilon}^*) dV = 2\pi e \int_{s_1}^{s_2} \left\{ r\Sigma_n \lambda^* + \left( \Sigma_t \cos \alpha + \Sigma_s \frac{r}{R_c} \right) u^* - \left( \Sigma_t \sin \alpha + \frac{d(r\Sigma_s)}{ds} \right) w^* \right\} ds + 2\pi e [w^* r \Sigma_s]_{s_1}^{s_2} \quad [7] $$

Internal pressure leads to the loading expression:

$$ \bar{F}_s = \Delta P \cdot \bar{n} \quad \text{on } \partial\Omega_f \quad [8] $$

Integration of external forces over $s_1$, $s_2$ gives the following term:

$$ \int_{0}^{2\pi s_2} \int_{s1} \Delta P \bar{n} \cdot \dot{V}^* r d\theta ds + \vec{F}_1 \cdot \vec{V}^* (s_2) + \vec{F}_2 \cdot \vec{V}^* (s_2) = 2\pi \Delta P \int_{s1}^{s2} \left( u^* - \frac{\lambda^* e}{2} \right) r ds + X_1 u_1^* + Z_1 w_1^* + X_2 u_2^* + Z_2 w_2^* \quad [9] $$

As well, the hydrostatic pressure term gives:

$$ \begin{aligned} \int_{\Omega} p \langle \mathrm{div} \, \vec{V}^* \rangle dV &= 2\pi e \int_{s1}^{s2} p \left( r\lambda^* + u^* \cos\alpha - w^* \sin\alpha + \frac{ru^*}{R_s} + r\frac{dw^*}{ds} \right) ds && [10] \\ &= 2\pi e \int_{s1}^{s2} \left( p r \lambda^* + p u^* \left( \cos\alpha + \frac{r}{R_s} \right) - w^* \left( p \sin\alpha + \frac{d(rp)}{ds} \right) \right) ds + 2\pi e [p w^*]_{s1}^{s2} \end{aligned} $$

Bringing these terms together, for any values of $\lambda^*$, $u^*$ and $w^*$ leads to local equations inside the membrane (i.e., for all “s”) and boundary conditions for $s_1$ and $s_2$:

$$ \begin{aligned} \Sigma_n - p &= -\frac{\Delta P}{2} & Z_1 &= -2\pi e (\Sigma_s - p)_{s_1} \\ (\Sigma_t - p) \cos \alpha + (\Sigma_s - p) \frac{r}{R_s} &= \frac{r\Delta P}{e} & Z_2 &= 2\pi e (\Sigma_s - p)_{s_2} \\ (\Sigma_t - p) \sin \alpha + \frac{d(r(\Sigma_s - p))}{ds} &= 0 & X_1 &= X_2 = 0 \end{aligned} \quad [11] $$

## Tube inflation : numerical and analytical comparison

The tube inflation problem is a particular case of the above problem. Since the membrane's profile (Fig.17) is a straight line $\sin\alpha$ vanishes as well as $1/R_s$ and $\cos\alpha$ equals to 1.

Local equations yields to the following system:

$$ \begin{aligned} \Sigma_n - p &= -\frac{\Delta P}{2} \\ (\Sigma_t - p) &= \frac{r\Delta P}{e} \\ (\Sigma_s - p) &= A \end{aligned} \quad [12] $$
Where A is a constant wrt. abscise s which is analogous to z in that case. The incompressibility condition gives:

$$ \lambda + \frac{u}{r} + w' = 0 \qquad [13] $$

<figure><img src="image_19.png" /><figcaption>Figure 17. Tube inflation problem</figcaption></figure>

The viscous behaviour is chosen to be Newtonian in this example:

$$ \begin{aligned} \Sigma_n &= 2\eta_0\lambda \\ \Sigma_t &= 2\eta_0 \frac{u}{r} \\ \Sigma_s &= 2\eta_0 w' \end{aligned} \qquad [14] $$

Since we have: u=dr/dt, w=dz/dt and λ=de/edt, we can mix the above relation in the following form:

$$ w' = - \left( \lambda + \frac{u}{r} \right) = - \left( \frac{\dot{e}}{e} + \frac{\dot{r}}{r} \right) = + \frac{\dot{L}}{L} \qquad [15] $$

$$ \begin{aligned} 2\eta_0 \frac{\dot{e}}{e} - p &= -\frac{\Delta P}{2} \\ 2\eta_0 \frac{\dot{r}}{r} - p &= \frac{r\Delta P}{e} \\ 2\eta_0 \frac{\dot{L}}{L} - p &= A \end{aligned} \qquad [16] $$

Introducing f as the aspect ratio (ie, f=e/r) we can solve the previous problem:

$$ f(t) + 2 = (f_0 + 2)\exp\left(-\frac{\Delta P}{4\eta_0}t\right) \qquad [17] $$

It worth noting that f(t) ratio is not influenced by the longitudinal elongation due to V. Since the length L(t) follows the relation [18]

$$ L(t) = L_0 + Vt \qquad [18] $$

The thickness and radius can be calculated from:
$$r(t) = \sqrt{\frac{e_0 r_0 L_0}{f(t) L(t)}} \qquad [19]$$

$$e(t) = f(t)r(t)$$

f(t) is a decreasing time function and the radius r(t) becomes infinite at ultimate time $t_U$ defined by $f(t_U) = 0$. Solving this relation gives:

$$t_U = \frac{4\eta_0}{\Delta P} \ln\left(\frac{f_0}{2} + 1\right) \approx \frac{2\eta_0 f_0}{\Delta P} \qquad [20]$$

The constant A can then be determined, pressure p and other stresses are then calculated:

$$\begin{align*} \sigma_n &= \sigma_r = -\frac{\Delta P}{2} \\ \sigma_t &= \sigma_\theta = \frac{r\Delta P}{e} \end{align*} \qquad [21]$$

$$\sigma_s = \sigma_z = 3\eta_0 \frac{\dot{L}}{L} + \Delta P \frac{2-f}{4f}$$

An adimensioned version of above equations can be obtained introducing reduced ratio f*, radius r*, thickness e* and length L*, as well as reduced time t*:

$$f^* = \frac{f(t)}{f_0}; R^* = \frac{R(t)}{R_0}; e^* = \frac{e(t)}{e_0}; L^* = \frac{L(t)}{L_0}; t^* = \frac{t(t)}{t_U} \qquad [22]$$

Using these definitions and introducing longitudinal strain rate from relation:

$$\dot{\varepsilon}_0 = \frac{V}{L_0} \Rightarrow \dot{\varepsilon}_0 = t_U \dot{\varepsilon}_0 \qquad [23]$$

The initial solution can be written as follow:

$$f^* = \frac{(f_0 + 2)}{f_0} \exp(-Ct^*) - \frac{2}{f_0} \quad \text{with:} \quad C = \ln\left(\frac{f_0 + 2}{2}\right) \qquad [24]$$

As a first comment of Fig.18, we will say that the assumption made on small values of aspect ratio is not compatible with $f_0$ values coming up to 1. The second comment is that the influence of $f_0$ is quite small on $f^*(t)$ evolution. This evolution is quasi linear from 1 to 0 à ultimate time $t_U$. Reduced radius and thickness are then given by:

$$r^* = \frac{1}{\sqrt{f(1 + \dot{\varepsilon}_0 t)}} \quad \text{and} \quad e = f^* r^* \qquad [25]$$

Fig.19 shows that elongation velocity have a significant influence on radius evolution. We can also say that initial aspect ratio has less and less influence since velocity increases. Figure 19 also shows thickness evolution which is the consequence of both previous plots. The same remark on aspect ratio $f_0$ can be made.

Solving equation [16] by a numerical technique, taking into account equation [15] leads to system:
$$p = -\eta_0 \frac{\dot{L}}{L} - \frac{\Delta P}{2} \left( \frac{r}{e} - \frac{1}{2} \right)$$

$$2\eta_0 \frac{\dot{e}}{e} = p - \frac{\Delta P}{2}$$

$$2\eta_0 \frac{\dot{r}}{r} = p + \frac{r\Delta P}{e}$$

[26]

<figure><img src="image_20.png" /><figcaption>**Figure 18.** *f*\* ratio evolution versus reduced time t\* for different initial aspect ratio f₀: 0, 0.1, ..., 0.9, 1</figcaption></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Reduced radius r\* and reduced thickness e\* evolutions versus reduced time t\* for different initial f₀ and different reduced strain rates: 0, 1, 10, 100</figcaption></figure>

Using a constant time step for iterative resolution, it appears that about 1000 time steps are needed for accurate resolution between initial time 0 and ultimate time t_u.
<figure><img src="image_23.png" /><figcaption>**Figure 20.** (a) Preform evolution, Newtonian case: uniform initial thickness $e_0 = 4$ mm and constant internal pressure $P_0 = 10^6$ Pa (b) Maximal radius and maximal length of the bottle versus time.</figcaption></figure>

Figure 20 shows the finite element simulation of free blowing using a Newtonian viscosity. An ultimate time $t_{U} = 0.09$ s can be highlighted. The preform shape becomes spherical and this is in contradiction with observed shapes of Fig. 16.

Accepted Manuscript