- Email: [email protected]
Large strain solutions of rubber components
~orn~u,~rs & Strucrurrs Pnnted I” Great Bntam
Vol
LARGE
17, No
5-6. PP 835-843.
STRAIN
@W-7949/83 Pergamon
1983
SOLUTIONS Bo
H~GBLAD
OF RUBBER
$3 00 + 00 Press Ltd.
COMPONENTS
and JAN ANDERS SUNDBERG
ASEA AB. S-721 83 VBsterHs, Sweden
Abstract-A
Mooney-Rlvlin
material model for plane stram and axisymmetric analyses of rubber has
been implemented in the ADINA computer code. Usmg a consistent penalty method, derived from a regularized mixed formulation. the nonhnear mcompressibility constraint IS weakened by a properly chosen projection procedure. The hydrostatic pressure (proportional to the penalty function and discontinuous at interelement boundaries) may assume constant or linear elementwlse variation, optimally matching the assumed displacement fields m four or nine node elements. respectively. Numerical studies on convergence rates and stability properties of pressure solutions give results that agree well with predictions made on the basis of the theory for the linear constraint case. The code has been used m the practical design of a highly strained rubber diaphragm in a sheet metal forming press. The pressure load was kept perpendicular to deformed surfaces and the growth of contact areas between rubber and metal boundaries durmg the forming process was simulated by introducing successive constramts in repeated restart runs.
I. INTRODUCTION Sheet metal components are often formed in fluid cell
presses where the metal is pressed against forming blocks by a thick rubber diaphragm which is backed up by oil pressure [l]. In each operating cycle the diaphragm is highly strained and it is necessary to optimize its geometrical shape in order to get an acceptable fatigue life. The ADINA-code [2] has been used in the analyses of the stress and strain state of the rubber diaphragm. The code contains as a standard option a two dimensional plane stress Mooney-Rivlin mode1 where it is possible to eliminate the pressure on element level and maintain exact incompressibility. However, for realistic analyses of the problem at hand it was necessary to implement a MooneyyRivlin model for axisymmetric and plane strain conditions. A penalty type formulation for the incompressibility constraint was chosen mainly because of its ease of implementation in an existing code and its supposed cost-effectiveness and good stability properties in comparison with mixed methods. Also the possiblity of simulating the compressive (near incompressible) behaviour of confined rubber components was of concern. Numerous investigators have contributed to the understanding of theoretical and numerical properties of the finite element penalty method applied to linear problems, see, e.g. [3-71. However, the precise conditions for robust numerical implementation of the methods are evidently not fully investigated [4. 51. For problems with nonlinear constraints, such as the incompressibility research on penalty
condition in finite elasticity. methods is in a preliminary
the
stage. The proposed methods are largely motivated by heuristic reasoning with the linear case as a model, see [8-121. In this work a generalization of one of the more promising penalty methods in the linear case, the so called consistent penalty method [5, 121 is applied to nonlinear elasticity. This method. where the con-
straint is weakened by a consistent projection procedure, may be viewed as a mixed method with a perturbed (regularized) Lagrangian. The projection procedure used in this work is essentially the same as outlined in Ref. [ 121.The penalty function (pressure) is projected elementwise onto a constant or linear polynomial function space. According to the linear theory this should permit optimal use of isoparametric elements with four or nine nodes, respectively. Some numerical studies have been performed on the stability properties and the convergence rates of the used finite elements. The convergence of the pressure in the &-norm boundary value problem
is investigated for a simple with a smooth solution field.
In analysing the engineering problem (the diaphragm-structure) a minor modification of the ADINA computer code was made to allow for pressure loading perpendicular to deformed surfaces. For comparison the problem was also analysed with a Hookean material model. using “effective” material constants. By introducing successive constraints in repeated restart runs it was possible to simulate the growth of contact areas between rubber and metal boundaries. 2. BASIC FORMULATION
We follow the exposition given in Ref. [ 131and use the same conventions of notation. The motion of the rubber body is studied in a stationary Cartesian coordinate system (x,, x2, .uJ. Using the so-called total Lagrangian formulation the principle of virtual displacements states that
s
I + bL,,Od v ’ + A;s,,6
=
r + A!R
(1)
0,
for the equilibrium in the configuration at time t + At. In (1) a left superscript indicates the configuration in which the quantity occurs. A left subscript indicates 835
836
Bo H~GBLAD
and JAN ANDERS SUNDBERG
the configuration with respect to which the quantity is measured. The integral in (1) is taken over the volume in the initial undeformed configuration at time t = 0. ’+*is,, is a Cartesian component of the 2nd PiolaKirchhoff stress tensor; ’ +*&, a Cartesian component of the Green-Lagrange strain tensor; and ‘+*‘R the virtual work of the externally applied loads and inertia forces. For hyperelastic materials the eqn (1) is highly nonlinear and an incremental Newton-Raphson procedure mav be needed for its numerical solution. Thus, we decompose in increments: ’ +A&s,, = &!?,,+ ,s,
3. A PENALTY
METHOD
We study incompressible, rubber-like characterized by a strain energy function
materials,
W(I,, 12)= C,(f, - 3) + C?(I? - 3) and a nonlinear incompressibility
(8)
constraint
G(I,) = 0
(9)
C, and C, are the Mooney-Rivlin constants. I,, I2 and of the Cauchy-Green deformation tensor, defined as
I, are the invariants
(2) ;c=1+2;c.
(10)
and ,+A,
Thus we have ocg = k,
(3)
+ atI/
I, = trace ;IC where
I, = i(trace AC) - i(trace bC*) bt,, = f (bu,, + ;a,,, + b%.,b&J)
(4a)
O5,= ()eg+ OCR
(4b)
A typical choice of the function G is
with
G(Z,) = I, - 1 oe,,= f(ou,, + or5, + &., ouA.,+ Ouk.r &.J
(5)
and ’ u!Y.,O~hJ or?,,= IO
(6)
bu,,ou,(i = 1,2,3) are the Cartesian components of the displacement field and their increments, respectively. Inserting (2)-(6) in (1) gives the equation for incremental quantities:
s
oS,~o~,Od V+
01
s
(11)
I3 = det $Z.
(12)
which inserted in (9) implies that all motions of the body are volume preserving. A classical way of handling a constraint like (9) is to change to a mixed formulation by introducing a Lagrange multiplier field ‘p that corresponds to the condition (9). The physical meaning of the multiplier will be a quantity that is proportional to the hydrostatic pressure. The displacement field $I and the pressure ‘p may be found as a stationary point (Au, ‘p) of the Lagrangian of the system. At this point the following variational relations are satisfied:
bS,,6ov,, “dV
OI
=
r+*‘R
_
s
s
&S,,G,e,OdV. (7)
01
0,
GFp’dV=O
(13)
This is the basic equation used in ADINA for total Lagrangian formulation.
‘R is the virtual work of external forces, at time t.
Remark
Regularize turbation
1 The virtual work ‘+*‘R contams contributions from pressure loads on external surfaces. As the pressure loads are defoimation dependent, ‘+ *‘R contains parts that may be transferred to the tangent stiffness matrix on the 1.h.s. of (7) (the so-called load stiffness matrix). Using equilibrium iterations (modified Newton-Raphson) it is reasonable to neglect this contribution and just update the load unbalance vector on the r.h.s. of (7). In general the load stiffness matrix is nonsymmetric. But in some cases it may be shown to be symmetric [14]. This indicates that the total pressure loading system is conservative. In this case, when using a symmetric quasi-Newton update of the tangent stiffness matrix (such as the BFGS-method in ADINA) we may hope for reasonable convergence even for relatively large pressure loads. For the actual problems discussed in this paper the pressure loading systems are indeed conservative.
--t
the system (13) by introducing
I
0,
‘pS’p OdV +
s 0)
GFp’dV=O
a per-
(15)
where c is a small quantity. In the continuous case (15) implies that the relation ‘p = fG(&)
(16)
is valid pointwise throughout the structure. In finite element analyses we choose the multiplier field ‘p to be discontinuous at interelement boundaries. Using (15) we can eliminate ‘p on the element level and reduce the total number of unknowns substantially. The relation (15) then shows that we are working with
837
Large strain solutions of rubber components a weakened constraint function Gh obtained by projecting G elementwise on.the space spanned by the chosen basis function for the pressure. Instead of (16) we get elementwise ‘p;I = ; Gh(Z3).
+o”T~[Gh(-&,~k~)~
(17)
When inserting (17) in (14) we also notice that we then are using a projected first variation of the constraint function
(23)
where
where Auare the nodal displacements of the elements and the subscript h denotes a projected quantity. Instead of (14) we get for an element e In (22a) and (23a) we have used the easily proved relations
ak Ok1 Ou-$'J= Ockl
(19) The eqn (19) is the analogue of (1) for an arbitrary element. It can be viewed as the condition for equilibrium of a body with a strain energy function defined elementwise as
This is realized by noticing that in this case the two operations “projecting” and “differentiating with respect to nodal degrees of freedom” commute. In (19) and (20) G,, is used as a penalty function with 6 as a penalty parameter. By chasing a penalty function with a reasonable physical behaviour we can also simulate near incompressibility in (20). In this work we use G(I,) = - f In (1J.
where
00
0
and
Inserting (22) and (23) in (19) we get the following Newton-Raphson scheme:
(21)
The logarithmic function conforms to measurements and has the right singular behaviour when Z, goes to zero. The stored energy tends to infinity if the body is compressed strongly. In finite elasticity, when assembling the element contributions in (19), we get a highly nonlinear system and some kind of Newton-Raphson procedure is usually needed in numerical computations. Denoting the increments of the nodal displacements and the Green-Lagrange strain tensor by ou and Ot,l respectively, we get in incrementing terms in (19):
+ZouT& G%J)
‘j-
(224
where the functions in the integrals on the left hand side are given by (22a) and (23a). Assembling the elementwise contributions in (26) we get a global relation involving a tangent stifiness matrix on the left hand side. Solving for the increments ou, the total displacements are updated according to ‘+~$l=&l+ou
(27)
after which (26) may be used in a modified Newton-Raphson scheme (equilibrium iterations in ADINA). The assembled form of (26) is the analogue of (7) in ADINA. In practical implementation it is convenient to adjust (26) so that zero pressure corresponds to zero strain. This is done by replacing G,,/E with G,,/c - 2C, -4C, in (26). Remark 2 The penalty method outlined above may be viewed as a straight forward generalization of the so-called consistent penalty method advocated in Ref. [5] for linear problems. For nonlinear elasticity various projecting methods are investigated in [g-12]. In the special case when the penalty function is projected on a constant space elementwise, we get the
Bo H~GGBLAD and
838
mean-dilatation formulation, already advocated in [8] for plasticity problems Remark
JAN ANDERS SUNDBERG
where M is the Gram-matrix a, =
3
It is well-known that in order to get useful numerical results with the penalty aproach for low order elements the penalty constraint has to be weakened. A popular weakening procedure is selective reduced integration [4,7]. In this work we use full integration and the weaking is accomplished and controlled by the choice of elemental basis functions for the pressure (the penalty function). Moreover, the used penalty method is by its construction equivalent to a mixed formulation.
r
Jh,
A4:dV
for &. 4, and & and (i = 0, 1, 2).
(31)
The integral of the product of two projected functions A, and B,, is computed elementwise according to
s
A,,&, OdV = (ao. aI, az) M-‘(& Ok?
B,. LUr
(32)
where 8, =
s
Bcj,‘dV
(i = 0, 1,2).
(33)
06,
Remark 4
For linear problems the usability of finite element penalty approach rests upon error estimates of the form (using simplified notations) [4,5]:
where I/ . /Iv and /I . lie denote norms on appropriate spaces for the displacement field (V) and the pressure (Q), respectively. d, and dz are constants, independent of E and the mesh size h. The values of the exponent m (typically 0.5 or 1) are determined by the regularity of the c-perturbation in (15). The exponent n is determined by the order of the assumed polynomial spaces for displacements and pressure. For a successful numerical implementation, the orders of the two spaces must be chosen so that a discrete Babuska-Brezzi condition is satisfied, see, e.g. Ref. [3,4]. For problems with nonlinear constraints the validity of an estimate like (28) is an open question. Numerical tests (see later on in the text) seems to indicate values around 1 for the exponent m when the solution field is smooth. We work with (28) as a hypothesis in our numerical experiments, performed to check convergence rates and stability properties of the used finite elements. 4. COMPUTATIONAL
ASPECTS
The notations in (22H27) agree with those used in ADINA [2, 131. From the structure of the formulas (22~(26) it is easily seen how the penalty function and its first and second order derivatives are combined with quantities already formed in the ADINA-code. The resulting 2nd Piola-Kirchhoff stress tensor is formed according to
where W is given by (8). The pressure is determined by (17). The projections in (23) and (23a) are efficiently performed at the same time as the “usual” unpenalized terms in (22) are integrated in ADINA. Suppose that the chosen pressure space is spanned by the basis function 40, C#J, and & (e.g. linear polynomial space). The projection of a function A is obtained as A,, = (bog & &)M-‘(aO, aI, a&T
(30)
In our implementation of the projecting procedure in ADINA the basis functions for the pressure can be specified on reference element level or globally. For the numerical tests, made in this paper, we have chosen to define the basis functions on the reference element level, where the degree of assumed polynomials for the displacements are unambiguously determined. It is believed that using globally defined basis functions for spanning, e.g. a linear pressure space may be non-optimal in a severely distorted nine-node isoparametric element. For regular meshes with rectangular elements the difference between the two methods of defining basis functions should vanish at least for small deformations. However, the question of the best way of specifying basis functions should be investigated further. An isoparametric element with four to nine nodes and constant or linear elementwise variation of pressure has been implemented in ADINA. According to the theory for linear constraints and uniform meshes [4] the following combinations are suitable for practical purposes (Table 1). Table
Element
I Pressure variation
4-node bilinear
constant
I-node serendipity
constant
9-node biquadratic
linear
The two last combinations also satisfy a discrete Babuska-Brezzi condition uniformly with decreasing mesh size h [4], and should therefore have good stability properties. In the nonlinear case the choice of the penalty parameter is a nontrivial task and may strongly influence the convergence of the solution procedure. For large strain analyses nonlinear iteration schemes, such as incremental Newton-Raphson, proceed by enforcing only linearisations of the constraint. The constraint is in general satisfied only by the converged solution at each increment. It is found that for large load increments an efficient way of reaching a solu-
Large
strain
solutions
tion is to imbed the problem in a series of problems with decreasing “compressibility” by using the penalty parameter as imbedding parameter. This is especially apparent in problems where the rubber is more or less confined by rigid boundary conditions. Then at the initial stage of an incremental process from an undeformed state, the tangent stiffness matrix only reflects the compressibility properties, i.e. the stiffness is dominated by l/c. Only after some iterations the material gradually begins to rearrange to allow for volume preserving deformations. Typically, when using ADINA, we start by incrementing up to full load with t(C, + C,) =O.Ol and then constrain c in two to ten steps per ten power decrease of t. The optimal final level of c is determined by a compromise between the actual accuracy of the model (discretization errors, modeling errors, etc.) and the available numerical precision of the computer. For an accurate model the discretization errors are small and hence the penalty parameter c should be chosen correspondingly small. In general smaller values of c(C, + C,) than 10e6 are seldom needed. 5. NUMERICAL
EXPERIMENTS
The influence of the numerical precision of the computer is illustrated in Fig. 1. Here the error in L,-norm between the penalty solution (pi) and the exact discrete solution (P,,) of a simple problem is plotted as a function of the penalty parameter. As seen from the figure the best accuracy is obtained when about half the numerical precision of the computer (8-9 decimal digits) is used for penalizing. In a practical situation the accuracy of the model is likely to motivate a much lower level of accuracy (see the horizontal line in the figure). If the Mooney-Rivlin constants deviate considerably from unity the quantity t(C, + C,) should be used instead of c along the horizontal axis in Fig. 1. The conclusions drawn from simple numerical tests such as illustrated in Fig. 1 are presumably quite general for numerical techniques using penalty methods. An indication of the regularity of t-perturbation in (I 5) can be obtained by studying the behaviour of the
of rubber
components
839
&-norm of the pressure solution as a function of 6 for some different kinds of problems. Figure 2 shows the behaviour of the &norm of the pressure when the penalty parameter E is decreasing for a rather highly strained case, defined in Fig. 2. The results in Fig. 2 suggest that for this smooth solution the exponent m in formula (15) has a value near 1. Similar experiments with a problem having a singular solution field (see Fig. 5) indicate a value of the exponent m that is much lower than 1. Figures 3 and 4 show numerical results from an analysis of the Chadwick-Haddon cylinder [ 151. The Mooney-Rivlin constants are: C, = 0.18134 N/mm’ C, = 0.0259 1 N/mm2 and the inner radius is doubled during the imposed deformation. The values at the Gauss-points are plotted. For the radial Cauchy stress component T,, good agreement with the theoretical results reported in [ 151 is obtained with relatively few elements. Both the stress and pressure solutions at Gausspoints of the 9-node element form a nearly smooth solution curve. This is also true for the mean-values taken at the centre of the 4-node elements. In order to test the stability properties of various isoparameteric elements when using the penalty method a problem with a singular solution field is analysed with 4-node bilinear, g-node serendipity and 9-node biquadratic elements, respectively. As has already been observed in the linear case [4], checkerboard pattern occurs in the use of 4-node elements. Notice also the stable appearance of the serendipity element (Fig. 5) with a constant pressure field. The pressure solution fields shown in Fig. 5 are far fom the solution of the incompressible case. The
stmn
Plane c,=15 C,=O5 Pressure IondIng -14’ L-node element lronstont pressure1
2
I
8
10
12
14
16
18
log ct1 occeptoble atrurocy level for o FE -model I” genera1
PIone
Max stretch ratlo F,nol L = 5 ~10~
Fig. I.
Dlsplarements are Imposed -1 3
Fig. 2.
stro1n
840
Bo H~GBLAD
and JAN ANDERS SUNDBERG
R [mm1 70,o
I
0 -
exact 4-node
X
n
9-node
qrHO~*
N/mm’1
I
I
I
I
I
98,O
,
” element element
(constant pressure) (Linear pressure)
’/
x
I
98 0
Fig.
pressure had not stabilized when c(C, + C,) was 10-4. The sensitivity to deviation from rectangular form has been studied by means of the axisymmetric meshes for the infinite cylinder shown in Fig. 6. The cylinder is subjected to internal pressure. The results from the test with 4-node elements with varying distortions are shown in the diagram. Of course increasing distortion means worse approximation
3.
properties but otherwise there is no tendency to oscillations in the pressure solutions. The 4-node element appears quite robust. Similar results are obtained for the 8- and 9-node elements (although the midside nodes or the center node were not moved in this case). Finally a simple convergence test has been performed to check the convergence rates of the pressure in 4- and 9-node elements (Fig. 7a). The same infinite
-4 lo14 -
element
,sL’
A
.05___-~_~-~~-=_~~__T--_ -10 -
12 -
-
A
P-node element
L-b-e
-
05 --
7 -L-node
-15 -
,z I
-
-----
_
----__
L-L-L__
pressure) e-node element (constantpressurel P-node element (hear pressure I 4 -node
element lconstont
4' /' +TI
10 lP' I /
a i
)R
980 Imml
70.0
Fig.
4.
L-node element
P-node element
Cl_15 C,.Ol66667
F,=OZ F,=O2
Fig.
5.
lpolnt IoBdl
__--
841
Large strain solutions of rubber components tz
L-node element
;
z
fconstont
pressure1
c,=1.5 c2=05 Pressure loodlng019367
R
t
:I
-P',
1
(lo'/ N/mm?
G
11
10
9-
7-
6,
LR 5
c, =1.5 c,=o.5 15
10
lmml
Pressure loading: 0.79367 Fig. 7(b).
Fig. 6.
In this case the finite element solution does not trivially satisfy the equilibrium conditions. The subdivision of the FE-mesh is indicated Fig. 7(b). Because of symmetry only a radial strip of quadratic elements is used in each
cylinder as in Fig. 6 is analysed.
refinement. The diagram confirms the theory that predicts h-convergence for 4-node elements and h2-convergence for 9-node in the &-norm for the pressure. 6. THE ENGINEERING PROBLEM
Sheet metal forming is a plastic forming process where a blank is pressed against a fixed die that determines the finished shape. The forming usually takes place without any significant thin-out of the blank. In fluid cell forming an oil volume, sealed by a rubber diaphragm, serves as a soft (universal) tool part, see Fig. 8. During the operation cycle the rubber diaphragm assumes the shapes of the rigid die and may therefore become highly strained. The finite element model of the diaphragmstructure is shown in Fig. 9. Four-node plane strain bilinear elements with elementwise constant pressure have been used. The growth of contact areas (see Fig. 9) between rubber and metal surfaces is simulated by
Fluid cell forming
o 9-node (Ilneor
010 h(d)
Fig. 7(a).
element pressure1
1
842
ESo H~GBLAD.
and JAN ANDERS SUNDBERG
deformed
Fig. 9.
7. CONCLUSIONS
introducing
successive constraints in repeated restart runs with ADINA. The results from different analyses in a crucial part of the structure are shown in Fig. 10. For comparison the problem was also analysed with a Mooney-Rivlin plane stress model and a Hookean material model with “effective” material constants. As seen from Fig. 10, the plane strain Mooney-Rivlin model gives results in rather good agreement with available results from measurements, while the plane stress Mooney-Rivlin model clearly is not applicable for this problem if correct quantitative results are wanted.
A Mooney-Rivlin material model for plane strain and axisymmetric analyses has been successfully implemented in ADINA using a consistent penalty projecting method. The numerical experiments indicate. that, as concerns the problems analysed in this work, this technique is generally applicable to nonlinear elasticity problems with smooth solution fields. For the non-smooth case (with point loads, re-entrant corners, etc.) the eight-node element with constant pressure and the nine-node element with linear pressure have a much more stable behaviour than the
stretch
ratio
undeformed
shape -.-
plane
strain
fiookean mtrl ” = 0.485 (reduced integration) plane Stress Mooney-Rivlln --_
Results
for
max.
stretch
ratio:
plane straul Hookean mtrl
1.48
plane strain Mooney-Rivl in
1.72
plane stress Mooney-Rivl in
1 .86
meas”rement
1.7-1.8
Fig. 10.
plane strain Mooney-Rivl in
Large
strain
solutions
element with constant pressure. These results agree well with what would be expected on the basis of the theory for the linear constraint case. On the other hand the simple four-node element is quite insensitive to distortions and its usefulness is demonstrated in the practical design of an engineering rubber component. The ADINA-code seems to be an ideal frame for implementing and testing new engineering facilities like those described in this work.
of rubber
four-node
Acknowledgements-The authors gratefully acknowledge the ASEA Company for the permission to publish this work and Prof. Klaus-Jiirgen Bathe, MIT, for valuable discussions and suggestions.
6.
7.
8.
9.
IO. REFERENCES
1. QUINTUS, sheet metal forming presses-a new means for metal forming with high pressures and single tools. ASEA J. 48, 141-146 (1975). 2. K. J. Bathe, ADINA: a finite element program for automatic dynamic incremental nonlinear analysis. Rep. 82448-1, Acoustics and Vibration Laboratory, Dept. of Mechanical engng, MIT. (1975). 3. M. Bercovier. Perturbation of mixed variational problems. Application to mixed finite element methods. R.A.I.R.O. (Numerical Analysis) 12, 21 l-236 (1978). 4. J. T. Oden and N. Kikuchi, Finite element methods for constrained problems in elasticity. Inf. J. Num. Meth. Engng. 18, 701-725 (1982). 5. M. S. Engelman. R. L. Sani, P. M. Gresho and M. Bercovter, Consistent vs reduced integration penalty methods for mcompressible media using several old and
CAS 'Id 17,No 5/&O
components
843
new elements. Int. J. Num. Meth. in Fluids 2. 2542 (1982). D. S. Malkus and T. J. R. Hughes, Mixed finite element methods-reduced and selective integration techniques: a unification of concepts. J. Comp. Meth. Appl. Mech. Engng. 15, 63-81 (1978). J. T. Oden, N. Kikuchi and Y. J. Song, Penalty-fimte element methods for the analysis of Stokesian flows. J. Camp. Meth. Appl. Mech. Engng. 31. 297-329 (1982). J. C. Nagtegaal. D. M. Parks and J. R. Rice, On Numerically accurate finite element solutions m the fully plastic range. J. Comp. Meth. Appl. Mech. Engng. 4, 153-157 (1974). T. J. R. Hughes, Generalization of selective integration procedures to anisotropic and nonlinear media. Inr. J. Num. Meth Engng 15, 1413-1418 (1980). D. S. Malkus. Finite elements with penalities in nonlinear elasticity. In!. J. Num. Mefh. Engng. 16. 121-136.
(1980). 11. E. Jankovich, F. Leblanc. M. Durand and M. Bercovier, A finite element for the analysis of rubber parts, expenmental and analytical assessment. Compuf. Structures 14, 385-391 (1981). 12. M. Bercovier, Y. Hasbani. Y. Gilon and K. J. Bathe, Finite element method analysis procedure for nonhnear incompressible elasticity. Proc. Inr. Symp. Hybrid and Mixed Finite Element Method, Atlanta (I 98 1). 13. K J. Bathe, Finite Element Procedures m Engmeering Analysis. Prentice-Hall. Englewood Cliffs, New Jersey (1982). 14. H. D. Hibbitt. Some follower forces and load stiffness. Int. J. Num. Meth. Enznz. 14. 937-941 (1979). of 15. P. Chadwick and E. WY Haddon, Inflation-extension a tube of incompressible isotropic elastic material. J. Inst. Maihs. Applies 10, 258-278 (1972).
Vol
LARGE
17, No
5-6. PP 835-843.
STRAIN
@W-7949/83 Pergamon
1983
SOLUTIONS Bo
H~GBLAD
OF RUBBER
$3 00 + 00 Press Ltd.
COMPONENTS
and JAN ANDERS SUNDBERG
ASEA AB. S-721 83 VBsterHs, Sweden
Abstract-A
Mooney-Rlvlin
material model for plane stram and axisymmetric analyses of rubber has
been implemented in the ADINA computer code. Usmg a consistent penalty method, derived from a regularized mixed formulation. the nonhnear mcompressibility constraint IS weakened by a properly chosen projection procedure. The hydrostatic pressure (proportional to the penalty function and discontinuous at interelement boundaries) may assume constant or linear elementwlse variation, optimally matching the assumed displacement fields m four or nine node elements. respectively. Numerical studies on convergence rates and stability properties of pressure solutions give results that agree well with predictions made on the basis of the theory for the linear constraint case. The code has been used m the practical design of a highly strained rubber diaphragm in a sheet metal forming press. The pressure load was kept perpendicular to deformed surfaces and the growth of contact areas between rubber and metal boundaries durmg the forming process was simulated by introducing successive constramts in repeated restart runs.
I. INTRODUCTION Sheet metal components are often formed in fluid cell
presses where the metal is pressed against forming blocks by a thick rubber diaphragm which is backed up by oil pressure [l]. In each operating cycle the diaphragm is highly strained and it is necessary to optimize its geometrical shape in order to get an acceptable fatigue life. The ADINA-code [2] has been used in the analyses of the stress and strain state of the rubber diaphragm. The code contains as a standard option a two dimensional plane stress Mooney-Rivlin mode1 where it is possible to eliminate the pressure on element level and maintain exact incompressibility. However, for realistic analyses of the problem at hand it was necessary to implement a MooneyyRivlin model for axisymmetric and plane strain conditions. A penalty type formulation for the incompressibility constraint was chosen mainly because of its ease of implementation in an existing code and its supposed cost-effectiveness and good stability properties in comparison with mixed methods. Also the possiblity of simulating the compressive (near incompressible) behaviour of confined rubber components was of concern. Numerous investigators have contributed to the understanding of theoretical and numerical properties of the finite element penalty method applied to linear problems, see, e.g. [3-71. However, the precise conditions for robust numerical implementation of the methods are evidently not fully investigated [4. 51. For problems with nonlinear constraints, such as the incompressibility research on penalty
condition in finite elasticity. methods is in a preliminary
the
stage. The proposed methods are largely motivated by heuristic reasoning with the linear case as a model, see [8-121. In this work a generalization of one of the more promising penalty methods in the linear case, the so called consistent penalty method [5, 121 is applied to nonlinear elasticity. This method. where the con-
straint is weakened by a consistent projection procedure, may be viewed as a mixed method with a perturbed (regularized) Lagrangian. The projection procedure used in this work is essentially the same as outlined in Ref. [ 121.The penalty function (pressure) is projected elementwise onto a constant or linear polynomial function space. According to the linear theory this should permit optimal use of isoparametric elements with four or nine nodes, respectively. Some numerical studies have been performed on the stability properties and the convergence rates of the used finite elements. The convergence of the pressure in the &-norm boundary value problem
is investigated for a simple with a smooth solution field.
In analysing the engineering problem (the diaphragm-structure) a minor modification of the ADINA computer code was made to allow for pressure loading perpendicular to deformed surfaces. For comparison the problem was also analysed with a Hookean material model. using “effective” material constants. By introducing successive constraints in repeated restart runs it was possible to simulate the growth of contact areas between rubber and metal boundaries. 2. BASIC FORMULATION
We follow the exposition given in Ref. [ 131and use the same conventions of notation. The motion of the rubber body is studied in a stationary Cartesian coordinate system (x,, x2, .uJ. Using the so-called total Lagrangian formulation the principle of virtual displacements states that
s
I + bL,,Od v ’ + A;s,,6
=
r + A!R
(1)
0,
for the equilibrium in the configuration at time t + At. In (1) a left superscript indicates the configuration in which the quantity occurs. A left subscript indicates 835
836
Bo H~GBLAD
and JAN ANDERS SUNDBERG
the configuration with respect to which the quantity is measured. The integral in (1) is taken over the volume in the initial undeformed configuration at time t = 0. ’+*is,, is a Cartesian component of the 2nd PiolaKirchhoff stress tensor; ’ +*&, a Cartesian component of the Green-Lagrange strain tensor; and ‘+*‘R the virtual work of the externally applied loads and inertia forces. For hyperelastic materials the eqn (1) is highly nonlinear and an incremental Newton-Raphson procedure mav be needed for its numerical solution. Thus, we decompose in increments: ’ +A&s,, = &!?,,+ ,s,
3. A PENALTY
METHOD
We study incompressible, rubber-like characterized by a strain energy function
materials,
W(I,, 12)= C,(f, - 3) + C?(I? - 3) and a nonlinear incompressibility
(8)
constraint
G(I,) = 0
(9)
C, and C, are the Mooney-Rivlin constants. I,, I2 and of the Cauchy-Green deformation tensor, defined as
I, are the invariants
(2) ;c=1+2;c.
(10)
and ,+A,
Thus we have ocg = k,
(3)
+ atI/
I, = trace ;IC where
I, = i(trace AC) - i(trace bC*) bt,, = f (bu,, + ;a,,, + b%.,b&J)
(4a)
O5,= ()eg+ OCR
(4b)
A typical choice of the function G is
with
G(Z,) = I, - 1 oe,,= f(ou,, + or5, + &., ouA.,+ Ouk.r &.J
(5)
and ’ u!Y.,O~hJ or?,,= IO
(6)
bu,,ou,(i = 1,2,3) are the Cartesian components of the displacement field and their increments, respectively. Inserting (2)-(6) in (1) gives the equation for incremental quantities:
s
oS,~o~,Od V+
01
s
(11)
I3 = det $Z.
(12)
which inserted in (9) implies that all motions of the body are volume preserving. A classical way of handling a constraint like (9) is to change to a mixed formulation by introducing a Lagrange multiplier field ‘p that corresponds to the condition (9). The physical meaning of the multiplier will be a quantity that is proportional to the hydrostatic pressure. The displacement field $I and the pressure ‘p may be found as a stationary point (Au, ‘p) of the Lagrangian of the system. At this point the following variational relations are satisfied:
bS,,6ov,, “dV
OI
=
r+*‘R
_
s
s
&S,,G,e,OdV. (7)
01
0,
GFp’dV=O
(13)
This is the basic equation used in ADINA for total Lagrangian formulation.
‘R is the virtual work of external forces, at time t.
Remark
Regularize turbation
1 The virtual work ‘+*‘R contams contributions from pressure loads on external surfaces. As the pressure loads are defoimation dependent, ‘+ *‘R contains parts that may be transferred to the tangent stiffness matrix on the 1.h.s. of (7) (the so-called load stiffness matrix). Using equilibrium iterations (modified Newton-Raphson) it is reasonable to neglect this contribution and just update the load unbalance vector on the r.h.s. of (7). In general the load stiffness matrix is nonsymmetric. But in some cases it may be shown to be symmetric [14]. This indicates that the total pressure loading system is conservative. In this case, when using a symmetric quasi-Newton update of the tangent stiffness matrix (such as the BFGS-method in ADINA) we may hope for reasonable convergence even for relatively large pressure loads. For the actual problems discussed in this paper the pressure loading systems are indeed conservative.
--t
the system (13) by introducing
I
0,
‘pS’p OdV +
s 0)
GFp’dV=O
a per-
(15)
where c is a small quantity. In the continuous case (15) implies that the relation ‘p = fG(&)
(16)
is valid pointwise throughout the structure. In finite element analyses we choose the multiplier field ‘p to be discontinuous at interelement boundaries. Using (15) we can eliminate ‘p on the element level and reduce the total number of unknowns substantially. The relation (15) then shows that we are working with
837
Large strain solutions of rubber components a weakened constraint function Gh obtained by projecting G elementwise on.the space spanned by the chosen basis function for the pressure. Instead of (16) we get elementwise ‘p;I = ; Gh(Z3).
+o”T~[Gh(-&,~k~)~
(17)
When inserting (17) in (14) we also notice that we then are using a projected first variation of the constraint function
(23)
where
where Auare the nodal displacements of the elements and the subscript h denotes a projected quantity. Instead of (14) we get for an element e In (22a) and (23a) we have used the easily proved relations
ak Ok1 Ou-$'J= Ockl
(19) The eqn (19) is the analogue of (1) for an arbitrary element. It can be viewed as the condition for equilibrium of a body with a strain energy function defined elementwise as
This is realized by noticing that in this case the two operations “projecting” and “differentiating with respect to nodal degrees of freedom” commute. In (19) and (20) G,, is used as a penalty function with 6 as a penalty parameter. By chasing a penalty function with a reasonable physical behaviour we can also simulate near incompressibility in (20). In this work we use G(I,) = - f In (1J.
where
00
0
and
Inserting (22) and (23) in (19) we get the following Newton-Raphson scheme:
(21)
The logarithmic function conforms to measurements and has the right singular behaviour when Z, goes to zero. The stored energy tends to infinity if the body is compressed strongly. In finite elasticity, when assembling the element contributions in (19), we get a highly nonlinear system and some kind of Newton-Raphson procedure is usually needed in numerical computations. Denoting the increments of the nodal displacements and the Green-Lagrange strain tensor by ou and Ot,l respectively, we get in incrementing terms in (19):
+ZouT& G%J)
‘j-
(224
where the functions in the integrals on the left hand side are given by (22a) and (23a). Assembling the elementwise contributions in (26) we get a global relation involving a tangent stifiness matrix on the left hand side. Solving for the increments ou, the total displacements are updated according to ‘+~$l=&l+ou
(27)
after which (26) may be used in a modified Newton-Raphson scheme (equilibrium iterations in ADINA). The assembled form of (26) is the analogue of (7) in ADINA. In practical implementation it is convenient to adjust (26) so that zero pressure corresponds to zero strain. This is done by replacing G,,/E with G,,/c - 2C, -4C, in (26). Remark 2 The penalty method outlined above may be viewed as a straight forward generalization of the so-called consistent penalty method advocated in Ref. [5] for linear problems. For nonlinear elasticity various projecting methods are investigated in [g-12]. In the special case when the penalty function is projected on a constant space elementwise, we get the
Bo H~GGBLAD and
838
mean-dilatation formulation, already advocated in [8] for plasticity problems Remark
JAN ANDERS SUNDBERG
where M is the Gram-matrix a, =
3
It is well-known that in order to get useful numerical results with the penalty aproach for low order elements the penalty constraint has to be weakened. A popular weakening procedure is selective reduced integration [4,7]. In this work we use full integration and the weaking is accomplished and controlled by the choice of elemental basis functions for the pressure (the penalty function). Moreover, the used penalty method is by its construction equivalent to a mixed formulation.
r
Jh,
A4:dV
for &. 4, and & and (i = 0, 1, 2).
(31)
The integral of the product of two projected functions A, and B,, is computed elementwise according to
s
A,,&, OdV = (ao. aI, az) M-‘(& Ok?
B,. LUr
(32)
where 8, =
s
Bcj,‘dV
(i = 0, 1,2).
(33)
06,
Remark 4
For linear problems the usability of finite element penalty approach rests upon error estimates of the form (using simplified notations) [4,5]:
where I/ . /Iv and /I . lie denote norms on appropriate spaces for the displacement field (V) and the pressure (Q), respectively. d, and dz are constants, independent of E and the mesh size h. The values of the exponent m (typically 0.5 or 1) are determined by the regularity of the c-perturbation in (15). The exponent n is determined by the order of the assumed polynomial spaces for displacements and pressure. For a successful numerical implementation, the orders of the two spaces must be chosen so that a discrete Babuska-Brezzi condition is satisfied, see, e.g. Ref. [3,4]. For problems with nonlinear constraints the validity of an estimate like (28) is an open question. Numerical tests (see later on in the text) seems to indicate values around 1 for the exponent m when the solution field is smooth. We work with (28) as a hypothesis in our numerical experiments, performed to check convergence rates and stability properties of the used finite elements. 4. COMPUTATIONAL
ASPECTS
The notations in (22H27) agree with those used in ADINA [2, 131. From the structure of the formulas (22~(26) it is easily seen how the penalty function and its first and second order derivatives are combined with quantities already formed in the ADINA-code. The resulting 2nd Piola-Kirchhoff stress tensor is formed according to
where W is given by (8). The pressure is determined by (17). The projections in (23) and (23a) are efficiently performed at the same time as the “usual” unpenalized terms in (22) are integrated in ADINA. Suppose that the chosen pressure space is spanned by the basis function 40, C#J, and & (e.g. linear polynomial space). The projection of a function A is obtained as A,, = (bog & &)M-‘(aO, aI, a&T
(30)
In our implementation of the projecting procedure in ADINA the basis functions for the pressure can be specified on reference element level or globally. For the numerical tests, made in this paper, we have chosen to define the basis functions on the reference element level, where the degree of assumed polynomials for the displacements are unambiguously determined. It is believed that using globally defined basis functions for spanning, e.g. a linear pressure space may be non-optimal in a severely distorted nine-node isoparametric element. For regular meshes with rectangular elements the difference between the two methods of defining basis functions should vanish at least for small deformations. However, the question of the best way of specifying basis functions should be investigated further. An isoparametric element with four to nine nodes and constant or linear elementwise variation of pressure has been implemented in ADINA. According to the theory for linear constraints and uniform meshes [4] the following combinations are suitable for practical purposes (Table 1). Table
Element
I Pressure variation
4-node bilinear
constant
I-node serendipity
constant
9-node biquadratic
linear
The two last combinations also satisfy a discrete Babuska-Brezzi condition uniformly with decreasing mesh size h [4], and should therefore have good stability properties. In the nonlinear case the choice of the penalty parameter is a nontrivial task and may strongly influence the convergence of the solution procedure. For large strain analyses nonlinear iteration schemes, such as incremental Newton-Raphson, proceed by enforcing only linearisations of the constraint. The constraint is in general satisfied only by the converged solution at each increment. It is found that for large load increments an efficient way of reaching a solu-
Large
strain
solutions
tion is to imbed the problem in a series of problems with decreasing “compressibility” by using the penalty parameter as imbedding parameter. This is especially apparent in problems where the rubber is more or less confined by rigid boundary conditions. Then at the initial stage of an incremental process from an undeformed state, the tangent stiffness matrix only reflects the compressibility properties, i.e. the stiffness is dominated by l/c. Only after some iterations the material gradually begins to rearrange to allow for volume preserving deformations. Typically, when using ADINA, we start by incrementing up to full load with t(C, + C,) =O.Ol and then constrain c in two to ten steps per ten power decrease of t. The optimal final level of c is determined by a compromise between the actual accuracy of the model (discretization errors, modeling errors, etc.) and the available numerical precision of the computer. For an accurate model the discretization errors are small and hence the penalty parameter c should be chosen correspondingly small. In general smaller values of c(C, + C,) than 10e6 are seldom needed. 5. NUMERICAL
EXPERIMENTS
The influence of the numerical precision of the computer is illustrated in Fig. 1. Here the error in L,-norm between the penalty solution (pi) and the exact discrete solution (P,,) of a simple problem is plotted as a function of the penalty parameter. As seen from the figure the best accuracy is obtained when about half the numerical precision of the computer (8-9 decimal digits) is used for penalizing. In a practical situation the accuracy of the model is likely to motivate a much lower level of accuracy (see the horizontal line in the figure). If the Mooney-Rivlin constants deviate considerably from unity the quantity t(C, + C,) should be used instead of c along the horizontal axis in Fig. 1. The conclusions drawn from simple numerical tests such as illustrated in Fig. 1 are presumably quite general for numerical techniques using penalty methods. An indication of the regularity of t-perturbation in (I 5) can be obtained by studying the behaviour of the
of rubber
components
839
&-norm of the pressure solution as a function of 6 for some different kinds of problems. Figure 2 shows the behaviour of the &norm of the pressure when the penalty parameter E is decreasing for a rather highly strained case, defined in Fig. 2. The results in Fig. 2 suggest that for this smooth solution the exponent m in formula (15) has a value near 1. Similar experiments with a problem having a singular solution field (see Fig. 5) indicate a value of the exponent m that is much lower than 1. Figures 3 and 4 show numerical results from an analysis of the Chadwick-Haddon cylinder [ 151. The Mooney-Rivlin constants are: C, = 0.18134 N/mm’ C, = 0.0259 1 N/mm2 and the inner radius is doubled during the imposed deformation. The values at the Gauss-points are plotted. For the radial Cauchy stress component T,, good agreement with the theoretical results reported in [ 151 is obtained with relatively few elements. Both the stress and pressure solutions at Gausspoints of the 9-node element form a nearly smooth solution curve. This is also true for the mean-values taken at the centre of the 4-node elements. In order to test the stability properties of various isoparameteric elements when using the penalty method a problem with a singular solution field is analysed with 4-node bilinear, g-node serendipity and 9-node biquadratic elements, respectively. As has already been observed in the linear case [4], checkerboard pattern occurs in the use of 4-node elements. Notice also the stable appearance of the serendipity element (Fig. 5) with a constant pressure field. The pressure solution fields shown in Fig. 5 are far fom the solution of the incompressible case. The
stmn
Plane c,=15 C,=O5 Pressure IondIng -14’ L-node element lronstont pressure1
2
I
8
10
12
14
16
18
log ct1 occeptoble atrurocy level for o FE -model I” genera1
PIone
Max stretch ratlo F,nol L = 5 ~10~
Fig. I.
Dlsplarements are Imposed -1 3
Fig. 2.
stro1n
840
Bo H~GBLAD
and JAN ANDERS SUNDBERG
R [mm1 70,o
I
0 -
exact 4-node
X
n
9-node
qrHO~*
N/mm’1
I
I
I
I
I
98,O
,
” element element
(constant pressure) (Linear pressure)
’/
x
I
98 0
Fig.
pressure had not stabilized when c(C, + C,) was 10-4. The sensitivity to deviation from rectangular form has been studied by means of the axisymmetric meshes for the infinite cylinder shown in Fig. 6. The cylinder is subjected to internal pressure. The results from the test with 4-node elements with varying distortions are shown in the diagram. Of course increasing distortion means worse approximation
3.
properties but otherwise there is no tendency to oscillations in the pressure solutions. The 4-node element appears quite robust. Similar results are obtained for the 8- and 9-node elements (although the midside nodes or the center node were not moved in this case). Finally a simple convergence test has been performed to check the convergence rates of the pressure in 4- and 9-node elements (Fig. 7a). The same infinite
-4 lo14 -
element
,sL’
A
.05___-~_~-~~-=_~~__T--_ -10 -
12 -
-
A
P-node element
L-b-e
-
05 --
7 -L-node
-15 -
,z I
-
-----
_
----__
L-L-L__
pressure) e-node element (constantpressurel P-node element (hear pressure I 4 -node
element lconstont
4' /' +TI
10 lP' I /
a i
)R
980 Imml
70.0
Fig.
4.
L-node element
P-node element
Cl_15 C,.Ol66667
F,=OZ F,=O2
Fig.
5.
lpolnt IoBdl
__--
841
Large strain solutions of rubber components tz
L-node element
;
z
fconstont
pressure1
c,=1.5 c2=05 Pressure loodlng019367
R
t
:I
-P',
1
(lo'/ N/mm?
G
11
10
9-
7-
6,
LR 5
c, =1.5 c,=o.5 15
10
lmml
Pressure loading: 0.79367 Fig. 7(b).
Fig. 6.
In this case the finite element solution does not trivially satisfy the equilibrium conditions. The subdivision of the FE-mesh is indicated Fig. 7(b). Because of symmetry only a radial strip of quadratic elements is used in each
cylinder as in Fig. 6 is analysed.
refinement. The diagram confirms the theory that predicts h-convergence for 4-node elements and h2-convergence for 9-node in the &-norm for the pressure. 6. THE ENGINEERING PROBLEM
Sheet metal forming is a plastic forming process where a blank is pressed against a fixed die that determines the finished shape. The forming usually takes place without any significant thin-out of the blank. In fluid cell forming an oil volume, sealed by a rubber diaphragm, serves as a soft (universal) tool part, see Fig. 8. During the operation cycle the rubber diaphragm assumes the shapes of the rigid die and may therefore become highly strained. The finite element model of the diaphragmstructure is shown in Fig. 9. Four-node plane strain bilinear elements with elementwise constant pressure have been used. The growth of contact areas (see Fig. 9) between rubber and metal surfaces is simulated by
Fluid cell forming
o 9-node (Ilneor
010 h(d)
Fig. 7(a).
element pressure1
1
842
ESo H~GBLAD.
and JAN ANDERS SUNDBERG
deformed
Fig. 9.
7. CONCLUSIONS
introducing
successive constraints in repeated restart runs with ADINA. The results from different analyses in a crucial part of the structure are shown in Fig. 10. For comparison the problem was also analysed with a Mooney-Rivlin plane stress model and a Hookean material model with “effective” material constants. As seen from Fig. 10, the plane strain Mooney-Rivlin model gives results in rather good agreement with available results from measurements, while the plane stress Mooney-Rivlin model clearly is not applicable for this problem if correct quantitative results are wanted.
A Mooney-Rivlin material model for plane strain and axisymmetric analyses has been successfully implemented in ADINA using a consistent penalty projecting method. The numerical experiments indicate. that, as concerns the problems analysed in this work, this technique is generally applicable to nonlinear elasticity problems with smooth solution fields. For the non-smooth case (with point loads, re-entrant corners, etc.) the eight-node element with constant pressure and the nine-node element with linear pressure have a much more stable behaviour than the
stretch
ratio
undeformed
shape -.-
plane
strain
fiookean mtrl ” = 0.485 (reduced integration) plane Stress Mooney-Rivlln --_
Results
for
max.
stretch
ratio:
plane straul Hookean mtrl
1.48
plane strain Mooney-Rivl in
1.72
plane stress Mooney-Rivl in
1 .86
meas”rement
1.7-1.8
Fig. 10.
plane strain Mooney-Rivl in
Large
strain
solutions
element with constant pressure. These results agree well with what would be expected on the basis of the theory for the linear constraint case. On the other hand the simple four-node element is quite insensitive to distortions and its usefulness is demonstrated in the practical design of an engineering rubber component. The ADINA-code seems to be an ideal frame for implementing and testing new engineering facilities like those described in this work.
of rubber
four-node
Acknowledgements-The authors gratefully acknowledge the ASEA Company for the permission to publish this work and Prof. Klaus-Jiirgen Bathe, MIT, for valuable discussions and suggestions.
6.
7.
8.
9.
IO. REFERENCES
1. QUINTUS, sheet metal forming presses-a new means for metal forming with high pressures and single tools. ASEA J. 48, 141-146 (1975). 2. K. J. Bathe, ADINA: a finite element program for automatic dynamic incremental nonlinear analysis. Rep. 82448-1, Acoustics and Vibration Laboratory, Dept. of Mechanical engng, MIT. (1975). 3. M. Bercovier. Perturbation of mixed variational problems. Application to mixed finite element methods. R.A.I.R.O. (Numerical Analysis) 12, 21 l-236 (1978). 4. J. T. Oden and N. Kikuchi, Finite element methods for constrained problems in elasticity. Inf. J. Num. Meth. Engng. 18, 701-725 (1982). 5. M. S. Engelman. R. L. Sani, P. M. Gresho and M. Bercovter, Consistent vs reduced integration penalty methods for mcompressible media using several old and
CAS 'Id 17,No 5/&O
components
843
new elements. Int. J. Num. Meth. in Fluids 2. 2542 (1982). D. S. Malkus and T. J. R. Hughes, Mixed finite element methods-reduced and selective integration techniques: a unification of concepts. J. Comp. Meth. Appl. Mech. Engng. 15, 63-81 (1978). J. T. Oden, N. Kikuchi and Y. J. Song, Penalty-fimte element methods for the analysis of Stokesian flows. J. Camp. Meth. Appl. Mech. Engng. 31. 297-329 (1982). J. C. Nagtegaal. D. M. Parks and J. R. Rice, On Numerically accurate finite element solutions m the fully plastic range. J. Comp. Meth. Appl. Mech. Engng. 4, 153-157 (1974). T. J. R. Hughes, Generalization of selective integration procedures to anisotropic and nonlinear media. Inr. J. Num. Meth Engng 15, 1413-1418 (1980). D. S. Malkus. Finite elements with penalities in nonlinear elasticity. In!. J. Num. Mefh. Engng. 16. 121-136.
(1980). 11. E. Jankovich, F. Leblanc. M. Durand and M. Bercovier, A finite element for the analysis of rubber parts, expenmental and analytical assessment. Compuf. Structures 14, 385-391 (1981). 12. M. Bercovier, Y. Hasbani. Y. Gilon and K. J. Bathe, Finite element method analysis procedure for nonhnear incompressible elasticity. Proc. Inr. Symp. Hybrid and Mixed Finite Element Method, Atlanta (I 98 1). 13. K J. Bathe, Finite Element Procedures m Engmeering Analysis. Prentice-Hall. Englewood Cliffs, New Jersey (1982). 14. H. D. Hibbitt. Some follower forces and load stiffness. Int. J. Num. Meth. Enznz. 14. 937-941 (1979). of 15. P. Chadwick and E. WY Haddon, Inflation-extension a tube of incompressible isotropic elastic material. J. Inst. Maihs. Applies 10, 258-278 (1972).