INTRODUCTION
The plastic industry often uses numerical simulation tools to reduce the time of the studies and to decrease the number of some expensive experimental tests. These tools are also used for optimising the process parameters of some known forming technologies as thermoforming and blow moulding of plastic bodies. The blowmoulding and/or thermoforming processes have the disadvantage, in the blowing step, to expand the hot parison until the contact with the cold mould, from where the sticking with this latter as the parison advances until this sticks completely to the mould. The consequence of this process is a product with an irregular distribution of thicknesses. In order to remedy this problem, the numerical simulation could be used for predicting the areas where the thicknesses are the thinnest and the value of the blowing pressure to apply.
Several plastic or rubber products are used in various modern technological branches. The modeling of the complex bahaviour of this kind of materials requires advanced mathematical and numerical tools (Mooney, 1940; Ogden, 1972; Crisfield, 1991). Nevertheless, to carry out a numerical simulation of the process such as thermoforming or blowing, we have to overcome several difficulties. The most awkward is the numerical incompressibility of the axisymmetric element and the large deformations of this kind of materials presumed hyperelastic or viscoelastic. The management of the contact between the parison and the mould (Heinstein et al., 2000; Zhong et al., 1994), with calculation of the reactions by a penalty method or by using Lagrange multipliers technique, is also a crucial difficulty to overcome.
A largely widespread approach for simulating this kind of materials is the MooneyRivlin’s model (Mooney, 1940; Rivlin, 1948), including as a special case the NéoHookean (Treloar, 1976) model where the function of strain energy density is expressed in terms of the invariants of the right CauchyGreen tensor. Another possibility is the use of the Ogden’s model, with a strain energy density formulated in terms of the principal elongations. Note that the Ogden model is considered as being one of the most complicated constitutive laws used in numerical modeling.
The main contribution of the present work comes from the observation that the
existing of finite element formulations for hyperelastic materials was not satisfactory
(Sze et al., 2004), particularly when using shell elements in large deformations.
Most of the efforts were concentrated on membrane elements which have the advantage
of having a plane stress behaviour. This allows eliminating incompressibility
problem in the finite element formulation. 3D volumetric elements (arbitrary
3D forms or 3D axisymetrical) (Ziane, 1999) have also numerical difficulties
due to the incompressibility and require a penalisation of the strain energy
density (Crisfield, 1991).
Recently, several works have been focused on the finite element developments. Sze developed a hyperelastic 18node solid element with a stabilization technique to eliminate both membrane and shear locking. Tuzel et al. (2004) studied the wave propagation in hyperelastic materials using Mooney Rivlin model for the axisymmetric problems. Düster et al. (2003) developed a higher order finite element (pversion) to treat hyperelastic materials, considered as quasiincompressible, which have been adapted for plane strain problems.
In this study, a new 3D axisymmetrical 4node element, labelled SFRQaxi (Space Fiber Rotation Quadrilateral Axisymmetric), has been developed, with a Mooney Rivlin hyperelastic model (Mooney, 1940), for axisymmetric problems. It’s based on the Space Fiber Rotation concept (Ayad, 1993, 2003) which allows introducing additional rotation degrees of freedom in order to give the element a high accuracy order without increasing the number of element nodes.
Two well known tests are presented. They show the good performances of the present element in terms of accuracy and computational time, when compared to the 8node quadratic element Q8axi for instance.
THEORETICAL FORMULATION OF THE ELEMENT SFRQAXI
The implementation of a hyperelastic law in the case of thin structures, with the assumption of plane stress hypothesis σ_{33} = 0, makes it possible to deal directly with the problem of incompressibility (I_{3 }= 1) and to determine the hydrostatic pressure p. In 3Dsolid or 3D axisymmetric formulations, the pressure cannot be determined by the boundary conditions. To overcome this difficulty, we adopt a quasiincompressible formulation suggested by Crisfield. It consists in uncoupling the strain energy in a volumetric energy, related to dilation and the energy of distortion, related to the isochoric deformation.
We present the main equations used for the implementation of the bahaviour law in an axisymmetric solid element. A Total Lagrangian Formulation, associated with a displacement model using the concept of Space Faber Rotation is used. The final model is a 4node 3D axisymmetric finite element with three degrees of freedom per node. 2 x 2 Gauss points are only used.
Approximation of the geometry: The position vector of an unspecified point in the meridian plan, expressed in a cylindrical base, is given by:
Formulation of the displacement model: The displacement field of an
unspecified point p(r, Z) of the plan is defined by the following approximations:
is the unit vector following the axis y perpendicular to plan (r, o, z), whereas
is
the position vector of the fiber connecting the point p to the node i (Fig.
1).
The term
(2) can be written under a matrix form:
Components ⟨u, w⟩ (2), can finally be written as:
The shape functions N_{i} are those of the classical 4node bilinear element Q4:
Derivatives of the terms (N_{i} (z–z_{i})) compared to R and Z are:
It should be noted that rotations θ_{y} are defined only with the nodes at the tops (Fig. 2).
Green Lagrange strain’s tensor: The strain gradient tensor components is written as:
The Right CauchyGreen tensor C = F^{T} F becomes

Fig. 1: 
Kinematics of a rotating fiber 

Fig. 2: 
Axisymmetric element SFRQaxi 
We define also the three invariants of:
The tensor of Green Lagrange is defined by:
Principle of the minimum total potential energy: The total potential energy is defined by:
with
W represents the strain energy per unit of no deformed volume, accumulated in the structure during the deformation:
Using a MooneyRivlin behaviour law, W given by Crisfield becomes:
Where
represent the modified invariants, whose expressions are given by (14):
The variation of the external potential energy is
given by:
⟨δu⟩ is the displacement in the local base
and {p} the compressive forces distributed along the meridian line.
Thus,
Where
represents the first variation of internal elementary energy:
Residual vector: The residual vector is defined from (17):
Using the MooneyRivlin law, we obtain
with
δ Ī_{1}, δ Ī_{2}, δ Ī_{3}
represent the variation of the invariants of tensor [C] (8):
Substituting the Eq. (21) and (17) leads
to following expression of the elementary internal strain potential:
with
By injecting the Eq. (23) in (22), becomes
with
and
The variation of the internal strain energy will be written as
with;
and
The vector of internal forces can be easily derived from
(25):
Tangent matrix (K^{T}): The calculation of the elementary tangent matrix can be explicitly defined by expressing the second variation of the elementary virtual work or by a perturbation technique. We have adopted the last technique which usually used. It is based on a finite difference scheme:
NUMERICAL EXPERIMENTS
All our geometry and results have been drawn and visualized using the educational
version of GID software (www.gidhome.com).
Circular plate under uniform pressure: We consider a circular plate treated by Hughes et al. (1983). Radius L and thickness h are respectively 7.5 in and 0.5 in. The plate is fixed on its circumference and subjected to a pressure on one of its faces (Fig. 3). We use a quasiincompressible hyperelastic material represented by the MooneyRivlin’s model with the constants C_{1} = 80 psi and C_{2} = 20 psi. The penalty factor K used to penalize the internal energy is taken equal to 10^{4}. Different meshes are used to analyse the performance of our SFRQaxi element in terms of CPU time, comparatively with the Q4axi element and the Q8axi element taken as reference solution.
Four different meshes are tested (Fig. 6). The last mesh, named same in the Fig. 6, allows getting the same number of degrees of freedom (dof) for the three elements: the problem is solved using the same memory.
Figure 5 shows that the results of the pressure variation, in accordance with the vertical displacement of the middle point of the plate, are almost identical compared to those given by Hugues et al. (1983). For a given strain step, we have represented on the Fig. 4 an image of the thickness distribution, showing in consequence the interest to use 3D axisymmetrical elements.

Fig. 3: 
Circular plate under pressure 

Fig. 4: 
Distribution of thicknesses 

Fig. 6: 
Element performancs in terms of CPU time 

Fig. 7: 
Hyperelastic circular membrane deformed by a punch 

Fig. 8: 
The membrane shape after deformation 

Fig. 9: 
Thickness distribution along radius 
Compared to the results of the quadratic element Q8axi in Fig. 6, our element SFRaxi performs well the Hughes solution with a small time. On the other hand, Q4axi seems turning more quickly. It cans give a no accurate solution indeed if the elements are distorted or the mesh is not sufficient.
Hyperelastic plate subjected to a rigid punch: A second example is a hyperelastic plate deformed by a rigid punch. This test has been experimentally studied by William (1970) for evaluating the thickness distribution of the plate (Fig. 7). A 2 x 20 finite element mesh has been used with the present SFRQaxi element and classical 4node Q4axi and 8node Q8axi elements. The numerical results are compared with those given by the experimental test (William, 1970) (Fig. 9).
Let’s consider the case of the sticking contact. We have used a contact algorithm based on the position code technique due to Heinstein et al. (2000). It will be detailed in a future issue. A simple NéoHookean model has been used to define the material behaviour. The rigidity modulus given by (William, 1970) is C_{1} = 0.448 MPa. The penalisation factor of the energy density is k = 100. The moving distance considered for the punch is 90 mm. An implicit scheme based on NewtonRaphson technique is used for the resolution. The test is carried out with 90 steps of vertical displacement (1 mm) applied to the punch (Fig. 8).
CONCLUSIONS
A formulation of a heuristic 3D axisymmetrical bilinear hyperelastic 4node element (12 dofs) for plastic forming processes is presented. For a similar accuracy, its performances compared with the expensive quadratic 8node Q8axi element (24 dofs) are better in terms of CPU time. The main advantage of the present element concerns the numerical integration of both tangent stiffness matrix and residual vector. Only 2x2 Gauss points have been used. Q8axi requires 3x3 Gauss points for instance. A second advantage is the possibility of carrying out remeshes at the areas of strong curvatures.
Moreover, compared to the membrane or shell elements which are defined by a middle neutral surface, SFRQ axi allows measuring directly the distribution of the thickness during the forming processes. An adequate penalisation term of the energy density is required in this case, in order to converge quickly and to give the anticipated results.
SFAQaxi can offer a good compromise between the execution speed and the required accuracy.