Abstract: This work deals with the non-linear formulation of 4-node axisymmetric hyperelastic solid model, based on the concept of SFR (Space Fiber Rotation). The SFRaxi model uses the kinematics of a space fibre to obtain an enriched displacement field. It improves in a significant way the precision of the classical 4-node axi- symmetrical element Q4axi. The corresponding results are comparable and even better in term of time CPU, with those of the higher order Q8axi element. A hyperelastic behaviour law based on the Mooney-Rivlin approach has been implemented allowing better simulations of blow-moulding and thermoforming of hollow plastic bodies. The numerical results, very promising, are mainly focused on some plastic forming tests without contact (swellings).
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 blow-moulding 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 Mooney-Rivlin’s model (Mooney, 1940; Rivlin, 1948), including as a special case the Néo-Hookean (Treloar, 1976) model where the function of strain energy density is expressed in terms of the invariants of the right Cauchy-Green 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 18-node 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 axi-symmetric problems. Düster et al. (2003) developed a higher order finite element (p-version) to treat hyperelastic materials, considered as quasi-incompressible, which have been adapted for plane strain problems.
In this study, a new 3D axi-symmetrical 4-node element, labelled SFRQaxi (Space Fiber Rotation Quadrilateral Axisymmetric), has been developed, with a Mooney Rivlin hyperelastic model (Mooney, 1940), for axi-symmetric 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 8-node 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 (I3 = 1) and to determine the hydrostatic pressure p. In 3D-solid or 3D axisymmetric formulations, the pressure cannot be determined by the boundary conditions. To overcome this difficulty, we adopt a quasi-incompressible 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 4-node 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:
| (1) |
Formulation of the displacement model: The displacement field of an unspecified point p(r, Z) of the plan is defined by the following approximations:
| (2) |
The term
| (3) |
Components ⟨u, w⟩ (2), can finally be written as:
| (4) |
The shape functions Ni are those of the classical 4-node bilinear element Q4:
| (5) |
Derivatives of the terms (Ni (z–zi)) compared to R and Z are:
| (6) |
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:
| (7) |
The Right Cauchy-Green tensor C = FT F becomes
| Fig. 1: | Kinematics of a rotating fiber |
| Fig. 2: | Axisymmetric element SFRQaxi |
| (8a) |
We define also the three invariants of:
| (8b) |
The tensor of Green Lagrange is defined by:
| (9) |
Principle of the minimum total potential energy: The total potential energy is defined by:
| (10) |
with
| (11) |
W represents the strain energy per unit of no deformed volume, accumulated in the structure during the deformation:
| (12) |
Using a Mooney-Rivlin behaviour law, W given by Crisfield becomes:
| (13) |
Where
| (14a-c) |
The variation of the external potential energy
| (15) |
⟨δu⟩ is the displacement in the local base
Thus,
| (16) |
Where
| (17) |
Residual vector: The residual vector is defined from
| (18) |
| (19) |
Using the Mooney-Rivlin law, we obtain
| (20) |
with
|
|
δ Ī1, δ Ī2, δ Ī3 represent the variation of the invariants of tensor [C] (8):
| (21) |
Substituting the Eq. (21) and (17) leads to following expression of the elementary internal strain potential:
| (22) |
with
| (23) |
By injecting the Eq. (23) in (22),
| (24) |
with
and
The variation of the internal strain energy will be written as
| (25) |
with;
| (26) |
and
| (27) |
| (28) |
The vector of internal forces can be easily derived from
| (29) |
Tangent matrix (KT): 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:
| (30) |
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 quasi-incompressible hyperelastic material represented by the Mooney-Rivlin’s model with the constants C1 = 80 psi and C2 = 20 psi. The penalty factor K used to penalize the internal energy is taken equal to 104. 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 axi-symmetrical elements.
| Fig. 3: | Circular plate under pressure |
| Fig. 4: | Distribution of thicknesses |
| Fig. 5: | Loading curve |
| 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 4-node Q4axi and 8-node 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éo-Hookean model has been used to define the material behaviour. The rigidity modulus given by (William, 1970) is C1 = 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 Newton-Raphson 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 axi-symmetrical bilinear hyperelastic 4-node element (12 dofs) for plastic forming processes is presented. For a similar accuracy, its performances compared with the expensive quadratic 8-node 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 re-meshes 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.