# Finite element analysis of non linear transversaly isotropic hyperelastic membranes for thermoforming applications

Nicolas Chevaugeon, Erwan Verron, Bernard Peseux

► To cite this version:

Nicolas Chevaugeon, Erwan Verron, Bernard Peseux. Finite element analysis of non linear transvers-
saly isotropic hyperelastic membranes for thermoforming applications. ECOMAS 2000, Oct 2000,
Barcelone, Spain. hal-01008224

HAL Id: hal-01008224

https://hal.science/hal-01008224

Submitted on 3 Nov 2016

**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
FINITE ELEMENT ANALYSIS OF NONLINEAR
TRANSVERSELY ISOTROPIC HYPERELASTIC
MEMBRANES FOR THERMOFORMING APPLICATIONS

N. Chevaugeon, E. Verron, and B. Peseux

Laboratoire Mécanique et Matériaux

Division Structures

École Centrale de Nantes

BP 92101

44321 Nantes cedex 3, France

e-mail: nicolas.chevaugeon@ec-nantes.fr

web page: http://www.ec-nantes.fr/Structures/

**Key words:** Hyperelasticity, Thermoforming, Fibers-reinforced material, Transversely isotropic constitutive equation

**Abstract.** Recent advances have been made to simulate plastic blow-moulding and thermoforming processes. However, most of the work done is focused on thermoplastic materials assumed to obey an isotropic behaviour. The present paper deals with free and confined inflation of heat-softened fibers-reinforced membranes in order to simulate moulding processes of composite thermoplastic parisons. The material is assumed to be hyperelastic and incompressible. To take into account the anisotropy induced by aligned fibers, the material symmetry class is considered transversely isotropic. A general expression of the strain energy function in terms of tensorial invariants is developed and governing equations of the inflation problem are solved by a total Lagrangian finite element method for space discretization and an explicit finite-difference algorithm for time integration. In order to illustrate our approach, both free and confined inflation examples are examined.
# 1 INTRODUCTION

The finite element simulation of thermoforming process for isotropic materials has been extensively examined in the past. Both axisymmetrical and three-dimensional models have been developed. In 1989, Zamani et al. [1] reviewed publications dealing with this subject. In the present introduction, we only focus on the more recent works on the subject. deLorenzi and Nied [2, 3] examine both blow-moulding and thermoforming processes simulation. Using experimental observations, they assume that the thermoplastic parison can be modeled by a rubberlike hyperelastic membrane. In the same period, Charrier and coworkers [4, 5] used the same approach for viscoelastic membranes. All these papers report on a quasi-static treatment of the problem. More recently, Bourgin et al. [6] and Verron et al. [7, 8] prefer to consider inertial effects and use a dynamic explicit numerical procedure to solve the problem.

All these works are focused on isotropic materials. There are just few studies concerning the thermoforming of fibers-reinforced thermoplastics. Some authors use the approach developed in injection molding process of fiber loaded material [9]. The material is described at a micromechanical scale (unit cell model) and is considered as a viscous fluid loaded by rigid fibers. From this microscopic point of view and using homogeneization technics, authors build macroscopic constitutive equations, formulated in terms of strain rate. The corresponding behaviour is anisotropic viscous fluid like [10, 11, 12].

The main goal of our work consists in using the classical developments relative to thermoforming of isotropic materials and in adapting them to the general case of anisotropic materials. Following the work of Spencer [13], we choose to describe the anisotropy induced by aligned fibers in the rubberlike thermoplastics at the macromechanical scale using material symmetry properties. The present paper reports on the use of transversely hyperelastic constitutive models in the thermoforming simulation. Such theory was successfully used by Kyriacou et al. [14] in membrane inflation problems.

In the next section the theoretical background used to describe stress and strain for finite deformation in the Lagrangian point of view is recalled, and the general form of the strain energy function for hyperelastic transversely isotropic materials is established. Section 3 is devoted to the presentation of the finite element formulation suitable to describe the inflation of membrane under an inflating pressure. The numerical implementation of the constitutive equation is highlighted. Some numerical results obtained by our method for free and confined inflation cases are presented in section 4. Last, concluding remarks and perspectives are proposed.

# 2 HYPERELASTIC TRANSVERSELY ISOTROPIC CONSTITUTIVE EQUATION

In this paper, fibers-reinforced thermoplastic materials are considered. Fibers are assumed to be aligned and materials are studied far above their glass transition temperature. Therefore materials obey hyperelastic transversely isotropic constitutive equations.
## 2.1 Stress and strain definitions

Consider a material point $P$ of a body $B$. $P$ occupies the position $X$ relative to a frame of reference $R_0 : \{O, e_i\}$ in the undeformed initial configuration $C_0$. In the deformed configuration $C$, $x$ is the position of $P$ relative to an observer frame of reference $R : \{o, e_i\}$. The transformation between $C_0$ and $C$ is defined by:

$$x = \xi(X, t) \qquad (1)$$

where $t$ is the present time.

In the neighborhood of $P$, the deformation gradient $F$, the right Cauchy-Green deformation tensor $C$ and the Green-Lagrange strain tensor $E$ are respectively defined by:

$$F(X, t) = \frac{\partial x(X, t)}{\partial X} \qquad (2)$$

$$C(X, t) = F^T F \qquad (3)$$

$$E(X, t) = \frac{1}{2} [C(X, t) - I] \qquad (4)$$

where $I$ is the second-order identity tensor and $\cdot^T$ denotes the transposition.

Let $\sigma$ and $S$ be the Cauchy "true" stress tensor defined in the deformed configuration and the second Piola-Kirchhoff stress tensor relative to the undeformed configuration. They are respectively defined by:

$$\mathbf{df} = \sigma \cdot \mathbf{n} \, ds \qquad (5)$$

$$\mathbf{dF} = S \cdot N \, dS \qquad (6)$$

in which $\mathbf{df}$ is the force in the deformed configuration acting on a deformed oriented area $\mathbf{n} \, ds$ and $\mathbf{dF}$ is the force in the undeformed configuration acting on an undeformed oriented area $N \, dS$. $\mathbf{df}$ and $\mathbf{n} \, ds$ are the transformed vectors of $\mathbf{dF}$ and $N \, dS$ by the deformation gradient.

With some standard manipulations the two stress tensors can be related by:

$$S = J F^{-1} \sigma F^{-T} \qquad (7)$$

where $J$ is the Jacobian of the transformation (for a volume-preserving deformation: $J = 1$).

## 2.2 Construction of a Strain energy function

Hyperelastic materials are defined by the existence of a scalar function $W(F)$ called stored energy or strain energy function, from which stresses can be derived at each point. In order to satisfy the objectivity principle, the strain energy function must be invariant under changes of observer frame of reference. It is well-known that the Cauchy-Green
deformation tensor is invariant under changes of observer frame of reference [15]. Thus, if the strain energy function W can be written as a function of C, it automatically satisfies the objectivity principle.

We now have to take into account the material symmetries that restrict the way W depends on C. Every orthogonal transformation member of the material symmetry group must keep the strain energy unchanged. For isotropic material, the material symmetry group consists in the entire group of orthogonal transformations. In the case of unidirectional fibers-reinforced composite, the symmetry group has to be determined. The fibers orientation is defined by a local unit vector denoted $\mathbf{a}_0(\mathbf{X})$. Therefore, the material symmetry group is the group of orthogonal transformations which keep $\mathbf{a}_0(\mathbf{X})$ unchanged. Consequently the material is isotropic in all planes orthogonal to the fibers direction and is called *transversely isotropic*. The general form of the corresponding strain energy function must satisfy:

$$W(\mathbf{C}) = W(\mathbf{Q} \cdot \mathbf{C} \cdot \mathbf{Q}^T) \quad \forall \mathbf{Q} \text{ orthogonal with } \mathbf{Q} \cdot \mathbf{a}_0(\mathbf{X}) = \mathbf{a}_0(\mathbf{X}) \quad (8)$$

Another way of writing W consists in introducing explicitly the dependence on $\mathbf{a}_0$ [13, 16]. The invariance with respect to the symmetry group is now expressed by:

$$W(\mathbf{C}, \mathbf{a}_0) = W(\mathbf{Q} \cdot \mathbf{C} \cdot \mathbf{Q}^T, \mathbf{Q} \cdot \mathbf{a}_0) \quad \forall \mathbf{Q} \text{ orthogonal} \quad (9)$$

As the sign of $\mathbf{a}_0$ has no physical meaning, we can introduce the orientation tensor $\mathbf{A} = \mathbf{a}_0 \otimes \mathbf{a}_0$ with no loss of information ($\otimes$ stands for the tensor outer product). Then the previous equation (9) becomes:

$$W(\mathbf{C}, \mathbf{A}) = W(\mathbf{Q} \cdot \mathbf{C} \cdot \mathbf{Q}^T, \mathbf{Q} \cdot \mathbf{A} \cdot \mathbf{Q}^T) \quad \forall \mathbf{Q} \text{ orthogonal} \quad (10)$$

Thus W is an isotropic function of two symmetric tensors C and A. Using representation theorems [17], W can be written as a function of the invariants of C, A and their products. Following equations give a possible set of ten independent scalar invariants on which W can depend [18]:

$$I_1 = \text{trC} \qquad I_2 = \frac{1}{2} [(\text{trC})^2 - \text{trC}^2] \qquad I_3 = \det C \quad (11)$$

$$A_1 = \text{trA} \qquad A_2 = \frac{1}{2} [(\text{trA})^2 - \text{trA}^2] \qquad A_3 = \det A \quad (12)$$

$$I_4 = \text{tr}(\mathbf{C} : \mathbf{A}) \qquad I_5 = \text{tr}(\mathbf{C}^2 : \mathbf{A}) \qquad I_6 = \text{tr}(\mathbf{C} : \mathbf{A}^2) \quad (13)$$

$$I_7 = \text{tr}(\mathbf{C}^2 : \mathbf{A}^2) \quad (14)$$

where $(I_i)_{i=1,3}$ are the invariants of C, $(A_i)_{i=1,3}$ those of A and $(I_i)_{i=4,7}$ are product invariants.
In the present case of transversely isotropic behaviour, $\mathbf{A} = \mathbf{A}^2$ and the invariants of $\mathbf{A}$ become:

$$A_1 = 1 \quad A_2 = 0 \quad A_3 = 0 \qquad (15)$$

Moreover, some relations can be written for the invariants of the products:

$$I_5 = I_7 \quad I_4 = I_6 \qquad (16)$$

Therefore, the number of relevant independent invariants on which $W$ depends is reduced to five [19]:

$$W = W(I_1, I_2, I_3, I_4, I_5) \qquad (17)$$

Note that in the isotropic case, the number of invariants is reduced to three [20]: $I_1, I_2$ and $I_3$.

Taking into account the incompressibility assumption, i.e.: $I_3 = 1$, the number of independent invariants is finally reduced to four.

The 2nd Piola-Kirchhoff stress tensor can be calculated by differentiation of the strain energy function with respect to the Cauchy-Green dilatation tensor $\mathbf{C}$ [15]:

$$\mathbf{S} = -p\mathbf{C}^{-1} + 2\frac{\partial W}{\partial \mathbf{C}} \qquad (18)$$

where $p$ is an arbitrary hydrostatic pressure preventing for any volume change. This pressure has no physical meaning and is determined from equilibrium equations. Using the differentiation chain rule, the stress tensor can be expressed as:

$$\mathbf{S} = -p\mathbf{C}^{-1} + 2 \sum_{i=1,2,4,5} \frac{\partial W}{\partial I_i} \frac{\partial I_i}{\partial \mathbf{C}} \qquad (19)$$

with:

$$\frac{\partial I_1}{\partial \mathbf{C}} = \mathbf{I} \qquad (20)$$

$$\frac{\partial I_2}{\partial \mathbf{C}} = I_1 \cdot \mathbf{I} - \mathbf{C} \qquad (21)$$

$$\frac{\partial I_3}{\partial \mathbf{C}} = I_3 \cdot \mathbf{C}^{-1} \qquad (22)$$

$$\frac{\partial I_4}{\partial \mathbf{C}} = \mathbf{a}_0 \otimes \mathbf{a}_0 \qquad (23)$$

$$\frac{\partial I_5}{\partial \mathbf{C}} = \mathbf{a}_0 \otimes \mathbf{C} \cdot \mathbf{a}_0 + \mathbf{a}_0 \cdot \mathbf{C} \otimes \mathbf{a}_0 \qquad (24)$$
# 3 FINITE ELEMENT FORMULATION

In this section, we briefly present the finite element formulation adopted to simulate the thermoforming process. For more details, reader can refer to [7, 8].

## 3.1 Governing equations

The polymeric sheet is supposed to be a membrane and is represented by a two-dimensional continuum. Its thickness becomes a function of the position on the mid-surface. Assuming that there are no body forces, the Principle of Virtual Work, is expressed on the undeformed configuration for the inertial effects and internal work, and on the deformed configuration for the pressure follower force work:

$$ \int_{B_0} \delta \mathbf{u} \cdot \rho_0 \ddot{\mathbf{u}} dV_0 = - \int_{B_0} \delta \mathbf{E} : \mathbf{S} dV_0 + \int_{\delta B} \delta \mathbf{u} \cdot \mathbf{T} dS \quad \forall \delta \mathbf{u} \quad (25) $$

where $B_0$ and $\delta B$ stand respectively for the volume of the undeformed membrane and the surface of the deformed membrane, $\rho_0$ is the mass density, $\ddot{\mathbf{u}}$ is the acceleration vector, $\mathbf{T}$ is the external surface force due to pressure, and $\delta \mathbf{u}$ is a virtual displacement vector compatible with displacement boundary conditions.

The contact with the mold is assumed sticky: the heated parison cools down and becomes stiffer after contact [3]. Thus the blowing pressure is not sufficient to deform the membrane. Every material points which enter in contact is fixed until the end of the simulation. This assumption simplifies the numerical procedure [21].

## 3.2 Spatial discretization

The previous equation (25) is discretized by the finite element method. In the membrane context, 3-nodes isoparametrical triangular elements are chosen. These elements only deform in their plane, remain flat and triangular. The use of these simple elements provides analytical integration which is time saving. But it has a major flow: lot of elements are needed to obtain a good stiffness in high curvature regions of the parison. To avoid this problem, an automatic remeshing technique described in [7] is implemented. After discretization and integration on each element surface, the problem reduces to a system of ordinary second order differential equations:

$$ [M]\{\ddot{U}(t)\} = F^{\text{ext}}(t) - F^{\text{int}}(t) \quad (26) $$

where $[M]$ stands for the mass matrix, $U(t)$ is the nodal displacement vector, $F^{\text{ext}}(t)$, and $F^{\text{int}}(t)$ are respectively the external and internal nodal force vectors. By the use of the special lumping technique [22], $[M]$ can be reduced to a diagonal matrix, splitting the system in $n$ independent differential equations where $n$ is the number of degrees of freedom of the model.
## 3.3 Time-integration scheme

In order to solve the previous system (26), we use the classical explicit second-order central difference method. The velocity of the degree of freedom $i$, $\dot{U}_i(t)$, and its acceleration, $\ddot{U}_i(t)$, are given by:

$$ \dot{U}_i(t) = \frac{U_i(t + \Delta t) - U_i(t - \Delta t)}{2\Delta t} \qquad (27) $$

$$ \ddot{U}_i(t) = \frac{U_i(t + \Delta t) - 2U_i(t) + U_i(t - \Delta t)}{\Delta t^2} \qquad (28) $$

where $\Delta t$ is the time step. Consequently, the problem (26) is reduced to the following system:

$$ U_i(t + \Delta t) = \frac{\Delta t^2}{M_{ii}^d} [F_i^{ext}(t) - F_i^{int}(t)] + 2U_i(t) - U_i(t - \Delta t) \qquad (29) $$

for each degree of freedom $i$.

It is well-known that the stability of this scheme depends on the time step $\Delta t$. $\Delta t$ must be less than the critical time step $\Delta t_{cr}$. $\Delta t_{cr}$ is the smallest time step required for a wave to cross an element. As the present system is highly non-linear (large strains and material behaviour), stiffness properties change during calculations. To ensure the convergence of the scheme, $\Delta t_{cr}$ is periodically estimated by determining the stiffest element in the whole mesh.

## 3.4 Constitutive equation implementation

In order to calculate the internal force term in Eq. (29), the relation between deformations and stresses must be computed in each element. Using the plane stress state assumption, deformation and stress matrices in a single triangular element have the following forms:

$$ [\mathbf{C}] = \begin{bmatrix} C_{11} & C_{12} & 0 \\ C_{12} & C_{22} & 0 \\ 0 & 0 & C_{33} \end{bmatrix} \qquad (30) $$

and:

$$ [\mathbf{S}] = \begin{bmatrix} S_{11} & S_{12} & 0 \\ S_{12} & S_{22} & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (31) $$

where 1, 2, 3 represent the local axis in the element plane. The direction 3 is the normal to the surface.
Using nodal initial positions and displacements in the global referential axis, $C_{11}$, $C_{12}$ and $C_{22}$ are easily computed in each element. From the incompressibility constraint $detC = 1$, we have:

$$C_{33} = \frac{1}{C_{11}C_{22} - C_{12}^2} \qquad (32)$$

By definition, $C_{33}$ is the square of the principal stretch ratio in the membrane thickness direction:

$$\sqrt{C_{33}} = \frac{h}{h_0} \qquad (33)$$

where $h$ and $h_0$ are respectively the deformed and undeformed thicknesses of the membrane. By now, all terms of the strain tensor can be computed in each element.

The components of $\mathbf{A}$ in element local axis are computed using the fibers direction $\mathbf{a}_0$ previously defined in each undeformed element:

$$[\mathbf{A}] = \begin{bmatrix} a_{01}^2 & a_{01}a_{02} & 0 \\ a_{01}a_{02} & a_{02}^2 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (34)$$

All terms corresponding to the third direction are null because the fibers always stay in the element plane.

For a given expression of the strain energy function, $\mathbf{S}$ can now be calculated in the element local frame. The hydrostatic pressure $p$ is determined by using:

$$S_{33} = 0 \qquad (35)$$

In order to illustrate this general method, a very simple strain energy function is chosen:

$$W = \Phi(I_1 - 3) + \Psi(I_2 - 3) + \Theta I_4 \qquad (36)$$

where $\Phi$, $\Psi$ and $\Theta$ are the material parameters. This constitutive equation is a transversely isotropic generalization of the widely used Mooney-Rivlin model [20]. Using Eqs (18), (32), (35) and (36), the hydrostatic pressure $p$ is given by:

$$p = 2 \frac{\Phi + \Psi(C_{11} + C_{22})}{C_{11}C_{22} - C_{12}^2} \qquad (37)$$

And the stress components becomes:

$$S_{11} = 2(-C_{22}C_{33}^2[\Phi + \Psi(C_{11} + C_{22})] + \Phi + \Psi(C_{22} + C_{33}) + \Theta A_{11}) \qquad (38)$$

$$S_{22} = 2(-C_{11}C_{33}^2[\Phi + \Psi(C_{11} + C_{22})] + \Phi + \Psi(C_{11} + C_{33}) + \Theta A_{22}) \qquad (39)$$

$$S_{12} = 2(C_{12}C_{33}^2[\Phi + \Psi(C_{11} + C_{22})] + \Phi + \Psi(-C_{12}) + \Theta A_{12}) \qquad (40)$$

where $C_{33}$ is given by Eq. (32).
# 4 NUMERICAL EXAMPLES

## 4.1 Validation: free inflation

In order to validate the implementation of the previous constitutive equation, our results are compared with those obtained by Kyriacou et al. [14]. They use a quasi-static approach with four nodes elements to solve membrane inflation problems.

The free inflation of an initially flat square membrane submitted to an internal pressure is examined. The length of the undeformed membrane is denoted $l_0$ and its thickness $h_0$. The inflating pressure is noted $P$. The material obey the modified Mooney-Rivlin model previously defined. To simplify the discussion the dimensionless inflating pressure $P^*$ is defined by:

$$P^* = \frac{Pl_0}{\Phi h_0} \qquad (41)$$

$\Phi$ is the first material parameter.

The square membrane is first biaxially pre-stretched in its undeformed plane until the deformed square length $l$ reaches $1.1 l_0$. Kyriacou et al. use this pre-stretching to avoid any singularity in the tangent stiffness matrix at the first loading step. It is to note that no pre-streching is needed to start the computation in our dynamical approach. Numerical values are set to:

$$\Phi = 1. \quad \Psi = 0. \quad \Theta = 1. \quad H = 10^{-3} \quad L = 1. \qquad (42)$$

In this example, the square membrane is meshed with 400 triangular elements. Kyriacou et al. use 100 four nodes finite elements. Moreover, when the fibers direction coincides with one of the global axis directions, i.e.: $\mathbf{a}_0 = \mathbf{X}$ or $\mathbf{a}_0 = \mathbf{Y}$, only a quarter of the membrane is meshed using symmetry properties.

First, the inflating pressure $P^*$ is increased from 1.0 to 9.5. Figure 1 shows the height of the center node of the square sheet versus the inflating pressure. At the maximum pressure, the height is equal to 0.52. Our results are identical to the previously published data. Figure 2 shows the height of the center node at $P^* = 6$. for different fibers directions. In the figure, a relative angle between one side of the square and the fibers direction equal to zero corresponds to $\mathbf{a}_0 = \mathbf{X}$. The curve has a sinus-like shape and is symmetrical with respect to the $45^\circ$ angle because of the geometrical symmetry of the problem. The membrane is less stiff when fibers are aligned with the diagonal direction of the square. The relative difference in height between our curve and Kyriacou's one does not exceed 0.03%.

It is well-known that the pressure-displacement curve corresponding to the free inflation of a rubberlike membrane exhibits limit points [23]. In order to be able to reach higher deformation ratios, the load is imposed using a constant gaz flow rate. The pressure at time $t$ is computed step by step by using the Perfect Gaz Law and the volume inside the
<figure><img src="image_2.png" /><figcaption>Figure 1: Height of the center node versus inflating pressure for a square membrane</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Height of the center node versus angle between the fibers direction and the square side for $P^* = 6$.</figcaption></figure>

membrane at $t$ [7]. Using this continuation technique, the limit point can be overcome as shown in Figure 3. Thus the membrane reaches higher deformations.

The deformed geometry of a transversely isotropic inflated square membrane for two different fibers directions is presented in Fig. 4. The lines represent the fibers direction. The classical geometrical symmetry of the isotropic case disappears: membranes tend to exhibit a bulge perpendicularly to the fibers direction.

Finally, this example demonstrates that our results show good agreement with those of Kyriacou *et al*. Moreover, by the use of a gaz flow rate loading instead of a pressure
<figure><img src="image_4.png" /><figcaption>Figure 3: Height of the center node versus inflating pressure for a square membrane loaded with a constant gaz flow rate</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: On top views, fibers are aligned with the diagonal of the square in the initial configuration. On bottom views, fibers are aligned with the X side of the square</figcaption></figure>

loading, our method is able to go beyond limit points.
## 4.2 Thermoforming application

In the present example, some results concerning a simple thermoforming application are examined. The inflation of an initial circular flat sheet in a cylindrical mold is considered. Here automatic remeshing is used in order to be able to fill the corners of the mold.

As mentioned elsewhere [3], one of the main goal of the numerical simulation of such processes is to predict the final thickness distribution of molded parts. In order to evaluate the influence of fibers-reinforcement on the final geometry, comparisons are made between the membrane thickness distributions in the isotropic and the transversely isotropic cases. Material parameters are set to:

$$ \Phi = 1. \quad \Psi = 0. \quad \Theta = 3. \qquad (43) $$

Figures 5 and 6 present respectively four inflation stage of the moulding of a cylinder in the isotropic and the anisotropic cases. Due to geometrical and material symmetries, only a quarter of the membrane is considered. In the isotropic case, the thickness distribution

<figure><img src="image_6.png" /><figcaption>Figure 5: Thermoforming stages of a circular plane sheet in a cylindrical box: isotropic case. Gray scale: relative membrane thickness</figcaption></figure>

is axially symmetric (see Fig. 5). Figure 6 highlights the loss of symmetry due to the presence of fibers. The fibers direction is represented by the lines on the parison. It is to note that in this case axisymmetrical calculations would be unable to predict the final thickness distribution.
<figure><img src="image_7.png" /><figcaption>Figure 6: Thermoforming stages of a circular plane sheet in a cylindrical box: transversely isotropic case. Gray scale: relative membrane thickness</figcaption></figure>
Figure 7 presents the membrane thickness evolution along two edges shown in the small figure in the corner. The results correspond to the previous simulations. In the

<figure><img src="image_8.png" /><figcaption>Figure 7: Membrane thickness distribution at the end of the process</figcaption></figure>

isotropic case, the same thickness distributions are observed along the lines XZ and YZ, and are represented by the black solid line in the figure. In the transversely isotropic case, fibers are oriented as in Figure 6 and distributions differ for the two lines. The results corresponding to lines XZ and YZ are respectively plotted with a pointed line and a dashed line. In both cases, isotropic and transversely isotropic (for the two lines), the highest thickness is observed near the clamped edge and the smaller thickness corresponds to the corner of the mold.

An important remark found in the bibliography [3] relative to the thermoforming of isotropic material concerns the membrane thickness distribution at the end of the process. It is said that this distribution does not depend on the material model. But as pointed out in the present example, the anisotropic material behaviour due to the presence of fibers highly influence thickness distribution in thermoformed parts.

# 5 CONCLUDING REMARKS

In this paper a general form of the strain energy function for transversely isotropic hyperelastic materials has been developed. This function depends on a limited set of tensor invariants and ensure the objectivity principle. Moreover as shown by Kyriacou et al., the formulation of the constitutive equation is easy implemented in a quasi-static finite element program. Here it is shown that the implementation in a dynamic explicit finite element code is easier. Efficiency of the method is demonstrated on both free and
confined inflation problems.

As a perspective, we are working on a viscoelastic counterpart of the present hypere-
lastic model and on models for misaligned thermoplastic composites using a more general
orientation tensor. In order to determine the material parameters, a biaxial rheometer
for high strain rates is in development.
REFERENCES

[1] N. G. Zamani, D. F. Watt and M. Esteghamatian, Status of the finite element method in the thermoforming process, *Int. J. Num. Meth. Eng.*, **28**, 2681-2693, (1989).

[2] H. G. deLorenzi and H. F. Nied, Blow molding and thermoforming of plastics: finite element modeling, *Comput. Struct.*, **26**, 197-206, (1987).

[3] H. G. deLorenzi and H. F. Nied, Finite element simulation of thermoforming and blow molding, In A. I. Isayev, editor, *Progress in Polymer Processing*, 117-171, Hanser Verlag, (1991).

[4] J. M. Charrier, S. Shrivastava and R. Wu, Free and constrained inflation of elastic membranes in relation to thermoforming - non-axisymmetric problems, *J. Strain Analysis*, **24**, 55-73, (1989).

[5] S. Shrivastava and J. Tang, Large deformation finite element analysis of non-linear viscoelastic membranes with reference to thermoforming, *J. Strain Analysis*, **28**, 31-51, (1993).

[6] P. Bourgin, I. Cormeau and T. Saint-Martin, A first step towards the modelling of the thermoforming of plastic sheets, *J. Mat. Processing Tech.*, **54**, 1-11, (1995).

[7] E. Verron, G. Markmann and B. Peseux, Dynamic inflation of non-linear elastic and viscoelastic rubberlike membranes, *Int. J. Numer. Meth. Engng*, submitted.

[8] G. Marckmann, E. Verron and B. Peseux, Finite element analysis of blow-moulding and thermoforming process using a dynamic explicit procedure, *Polym. Engng Sci.*, submitted.

[9] D. Bhattacharyya et al., *Composite sheet forming*, Composite Materials Series, **11**, Elsevier, (1997).

[10] R.K. Okine, A. J. Beaussart and R. B. Pipes, Numerical modeling of sheet forming processes for thermoplastic composites, *Composite Mat. Tech.*, **III**, 3-19, (1991).

[11] P. De Luca, A. K. Pickett, T. Queckborner and E. Haug, An explicit finite element solution for the forming prediction of continuous fibre-reinforced thermoplastic sheets *Composites Manufacturing*, **6**, 237-243, (1995).

[12] A. K. Pickett and A. F. Johnson, Numerical simulation of the forming process in long fibre reinforced thermoplastics, *Composite Mat. Tech.*, **V**, 233-242, (1996)

[13] A.J.M. Spencer, *Continuum Theory of the mechanics of fibre reinforced Composite*, Springer-Verlag, (1984)
[14] S. K. Kyriacou, C. Schwab and J. D. Humphrey, Finite element analysis of nonlinear orthotropic hyperelastic membranes, *Comput. Mech.*, **18**, 269-278, (1996).

[15] M. F. Beatty, Topics in finite elasticity: hyperelasticity of rubber, elastomers and biological tissues - with examples, *Appl. Mech. Rev.*, **40**, 1699-1734, (1987).

[16] S. Govindjee, J. A. Weiss and B. N. Maker, Finite implementation of incompressible, transversely isotropic hyperelasticity, *Comput. Methods in Appl Mech. Engng*, **135**, 107-128, (1996).

[17] J. M. Ball, Convexity conditions and existence theorems in nonlinear elasticity, *Arch. Rat. Mech. Anal.*, **63**, 337-403, (1977).

[18] S. Sidoroff, *Cours sur les grandes transformations*, Technical Report 51, Greco, (1982).

[19] A. J. Burton and J. Bonet, A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations, *Comput. Methods in Appl. Mech. Engng*, **162**, 151-164, (1998).

[20] L. R. G Treloar, The mechanics of rubber elasticity, *Proc. R. Soc. Lond.*, **A351**, 301-330, (1976).

[21] R. E. Khayat, A. Derdouri and A. Garcia-Réjon, Multiple contact and axisymmetric inflation of a hyperelastic cylindrical membrane, *J. Mech. Eng. Sci.*, **207**, 175-183, (1993).

[22] R. L. Taylor and O. C. Zienkiewicz, *The finite element method. Volume I: Basic formulation and linear problems*, MacGraw-Hill Compagny, fourth edition, (1994).

[23] R. E. Khayat, A. Derdouri and A. Garcia-Réjon, Inflation of an elastic cylindrical membrane: non-linear deformation and instability, *Int. J. Solids Structures*, **29**, 69-87, (1992).