- Email: [email protected]
A phenomenological constitutive model for rubberlike materials and its numerical applications
International Journal of Plasticity 19 (2003) 1019–1036 www.elsevier.com/locate/ijplas
A phenomenological constitutive model for rubberlike materials and its numerical applications D. Besdo, J. Ihlemann* Institut fu¨r Mechanik, Universita¨t Hannover, Appelstr. 11, 30167 Hannover, Germany
Abstract A new phenomenological inelastic constitutive model for rubberlike materials is presented and its good correspondence to cyclic measurements is demonstrated for both uniaxial tension- and simple shear-tests up to large deformations. The model implies strong nonlinearities, hysteresis, the influence of the loading history and remaining deformations after unloading. Results of finite element calculations are represented to show the suitability of the constitutive model within this method for practical applications. # 2003 Elsevier Science Ltd. All rights reserved. Keywords: Rubber; Constitutive model; Inelasticity; FEM; Large deformation
1. Introduction Numerous technical applications of rubberlike material lead to large deformations. Furthermore, inelastic effects are often of great significance to the function of rubber containing structures. To simulate the mechanical behavior of those continua usually the finite element method is applied, in which by now the purely numerical problems due to rubber properties like strong nonlinearity or vanishing compressibility are controllable fairly well. On the other hand there is still a need for suitable three-dimensional constitutive models. Up to now the molecular foundations of the inelastic behavior of filled vulcanizates under large deformations, so far as known, are not yet transferable to suitable constitutive equations for the FEM in a satisfactory manner.
* Corresponding author. Tel.: +49-511-762-4120; fax: +49-511-762-4182. E-mail address: [email protected] (J. Ihlemann). 0749-6419/03/$ - see front matter # 2003 Elsevier Science Ltd. All rights reserved. PII: S0749-6419(02)00091-8
1020
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
Thus, phenomenological constitutive models are the remaining tools to make simulations possible. To consider the material behavior phenomenologically, suitable and extensive measurements are needed to identify those characteristics which are to be modeled (as an example see Khan and Lopez-Pamies, 2002). Preferably, the boundary conditions of those experiments should be as easy and unique as possible. Homogeneous deformations are preferred to allow the precise consideration of the material behavior itself instead of geometrical influences. The constitutive model, which will be presented in this context, is intended to be useful to the simulation of loading processes with rubber–typical large deformations. Moreover, strong nonlinearities, hysteresis and softening effects caused by the loading history are to be included in the model to allow the consideration of stress distributions within complex structures at the end. For that purpose, special cyclic tension and shear tests with varying upper elongation bounds have turned out to be useful. In Fig. 1 the results of an uniaxial tension test and a symmetrical simple shear test are outlined, both carried out by the Dutch company Akzo Nobel with a carbon black filled rubber used in air springs. Starting each with an initially unloaded specimen, both experiments consist of five increasing upper elongation bounds with several cycles at each level, approaching stationarity. As desired, the interesting characteristics of the material behavior are easy to see.
2. Constitutive model Constitutive descriptions are expressed either in Eulerian or in Lagrangean denotation. Unfortunately, the first Piola–Kirchhoff-stress tensor T, whose coefficients are used in Fig. 1, does not belong to one of these denotations. Hence, a transformation into Cauchy-stresses or second Piola–Kirchhoff-stresses T~ is necessary. For the corresponding transformation of the experimental data, an ideal incompressibility is assumed. The stationary cycles of the tension-test in terms of second Piola–Kirchhoff-stresses are plotted in Fig. 2. Obviously, the curvature of the lines is accentuated. The loading as well as the unloading curve have a changing sign of the curvature. In an Eulerian description (Fig. 3) using Cauchy-stresses no change of the curvature sign can be recognized. Therefore, this description seems to be simpler and, hence, to be preferred for the development of a constitutive model. In the case of the shear-test the curves for the Eulerian description do not differ from the corresponding graph in Fig. 1. However, the final aim will be a formulation using the quantities of the Lagrangean denotation, since most finite element codes use a Lagrangean formulation of the constitutive equations for implementation. Moreover, in order to implement the constitutive model, of course it must be formulated in a general tensorial 3D-scheme. A general problem of simulating rubber behavior is the near-incompressibility of those materials. Typically, the ratio of the bulk modulus to the shear modulus is several thousand. Because of this, the material behavior concerning pure changes in volume and on the other hand concerning deformations without any density
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1021
Fig. 1. Results of cyclic uniaxial tension and symmetric simple shear tests carried out by the Dutch company Akzo Nobel (transient cycles are plotted only partially). First Piola–Kirchhoff-stress tensor T vs. deformation gradient tensor F.
1022
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
changes preferably should be modeled and calculated separately. The volume ratio and thus, the ratio of the reference density to the current density, is connected with J3, the third main invariant of the deformation gradient F, which is equal to the determinant of the matrix of the coefficients of F in any orthonormalized coordinates. J3 ¼ I3 ðFÞ ¼ det½Fab With the help of J3 modifications of usual quantities are calculable, which are insensitive to the magnitude of the volume ratio (Sussmann and Bathe, 1987). Those modified quantities will be characterized by a superscript M. The use of these modified quantities corresponds to the idea of the ‘surface of natural density’ described by Makosey and Rajagopal (2001). In this context, corresponding modifications of
the right Cauchy–Green tensor C and the Jaumann rate b of the left Cauchy–Green tensor b are needed (0 denotes the purely mathematical deviatoric part). M
C¼
J32=3 C;
0 M 1 b ¼ b b b;
M
M
b ¼ J32=3 b
The new constitutive model should be able to describe the stationary cycles and the transition to higher strain levels as precisely as possible without a too high number of material constants. Nevertheless, the hysteresis, the remaining strains
Fig. 2. Stationary cycles of the tension test in Lagrangean description. Second Piola–Kirchhoff-stress tensor T~ vs. Green–Lagrange strain tensor .
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1023
after unloading and the influence of the strain history should be modeled in satisfactory approximation. The strain history is fitted by a scalar parameter. Thus, the influence of old loading directions leading to a material anisotropy is completely neglected. To define the one parameter, a special invariant is used, which represents the maximum difference of two eigenvalues. XT ¼ maxfjXI XII j; jXII XIII j; jXIII XI jg Applied to the Cauchy-stresses this is the reduced stress corresponding to Tresca’s yield criterion. The scalar parameter representing the strain history, is to be the maximum value of this invariant of the volume ratio insensitive modifications of the Cauchy–Green tensors in the past. M M H CT ðtÞ ¼ max CT ð Þ; 0 4 4 t M
In the following CH T is called the history function. Several constitutive models for rubber have been developed using rheological models with serial and parallel combinations of elastic, viscous and plastic elements
Fig. 3. Stationary cycles of the tension test in Eulerian description. Cauchy stress tensor vs. left Cauchy–Green tensor b.
1024
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
(Bergstro¨m and Boyce, 1998; Kaliske and Rothert, 1999; Lion, 2000; Lion and Haupt, 2001; Lubarda et al., 2002; Lulei and Miehe, 2001; Reese and Govindjee, 1998). Instead of this, the formulation of the present paper is in general based on the fundamental idea, to approximate the measured cycles by two limiting lines as asymptotic curves. The changes of the stresses in the range prescribed by these limits, which will actually be limiting tensors in the model, may be described by a specific differential equation, while strictly speaking the complete stresses will be an additive composition of some basic stresses and the solution of this differential equation (see Fig. 4). Additional experiments with varying strain rates showed that although the succession of the deformations has a strong influence on the material behavior, the influence of the amount of the velocity is very small over a wide velocity range. Therefore, the material description is designed to be independent of the rate in the above mentioned sense. Due to this feature this constitutive model is called an intrinsic model. Those rate independent models are often used to describe the behavior of metals undergoing finite deformations (Houlsby and Puzrin, 2000; Makosey and Rajagopal, 2001; Puzrin and Houlsby, 2001), whereas many other materials, e.g. teflon (Khan and Zhang, 2001), exhibit a distinctly time dependent behavior.
Fig. 4. Computed stationary tension cycle with corresponding basic and limiting stresses.
1025
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
The stresses coming from the solution of the differential equation, providing the determinant part of the hysteresis, are identified by means of the (additional) second Piola–Kirchhoff-stresses T~ A . The remaining parts of the total stresses are the so called basic stresses, which are similar to neo Hookean stresses but not identical with them, because the parameter is not constant in this case, but it is a definite function of three material constants p1, p2, p3 and in particular the history function M
CH T ðtÞ. Altogether, the stresses are divided into one part, which is physically corresponding to the deviatoric part of the Cauchy stress tensor and another part, representing the hydrostatic material behavior containing the bulk modulus K. 0 M A ~ ~ T ¼ 2C þ T C C1 þ K J3 ðJ3 1ÞC1 with :
p2 ¼ p1 þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi M
2
1 þ p3 CH T
Supposing a totally incompressible material, the bulk modulus tends to infinity and the value of the volume ratio J3 is identical with one in any case. Thus, the hydrostatic pressure is indeterminate and has to be calculated directly from the equilibrium equations. The differential equation for the additional stresses T~ A should cause these stresses A ~ T to tend to the current limiting stresses T~ L . This requirement leads to an
expression for T~ A , the simple Lagrangean rate of T~ A . T~ A ¼ LT T~ L T~ A
M
is again a definite function of material constants and the history function CH T ðtÞ (see below). LT is an invariant of a Lagrangean equivalent to an Eulerian time derivative preserving the required intrinsic character of the model. The tensor L will
be introduced below. The combination of the time derivative T~ A and the time dependent invariant LT cause a similar effect as the use of the ‘inelastic strain pathlength’ of Makosey and Rajagopal (2001). Beyond the frame of the intrinsic description a relaxation term can be easily inserted if required. T~ A ¼ LT T~ L T~ A p9 T~ A
For the final version of the differential equation, an additional term is needed to prevent physically hydrostatic parts of the limiting stresses influencing the evolution of the physically deviatoric parts of the stresses T~ A .
1026
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036 1 T~ A ¼ LT T~ L T~ A p9 T~ A þ T~ A C C1 3 p4 with : ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi M 1 þ p3 CH T
As mentioned above, the formulation of the limiting stresses will be easier in the Eulerian description. These stresses are explained as a tensorial function (Truesdell and Noll, 1965) of the Jaumann rate of the left Cauchy–Green tensor (for a discussion of useful objective corotational rates see Xiao et al., 2000, 2001). To retain the intrinsic character of the description, this time derivative is divided by its own M
invariant b T . Moreover, the vicinity to the prestraining maximum is included with M
the help of the ratio of the current value of b T with its maximum value in the past, which is equal to the history function. 2 L ¼
00
M M
11
M
3
BB CC 7 16 6 expBBp7 b b T CC þ p8 b 7 @ @ A A 4 5 M M M J3 bH T bT bT
M
with :
M
b T ¼ CT
To transfer Mthis equation identically into the Lagrangean notation, a Lagrangean M
equivalent to b is introduced, which has the same eigenvalues as b itself. However, it should be mentioned, that this tensor L may be unsymmetric, but it is still suitable to be an argument of a tensorial function. 0 M 1 1 1 C C þ C C L¼ C 2 With the help of L the transfer to a pure Lagrangean notation is possible without further problems. Once more, is a function of material constants and the history M
function CH T ðtÞ.
2
3 11 M L C L TA A þ p8 5 C1 T~ L ¼ 4 exp@@p7 LT M H LT CT 0 with :
00
1
B C M B C p6 B1 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C ¼ p5 C H TB 2 C @ A M p26 þ CH T
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
Fig. 5. Tension and shear test: comparison of experimental and calculated data.
1027
1028
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
With the explanation of the limiting stresses the constitutive model is completed and ready for material behavior simulation and implementation in the finite element method. Apart from the bulk modulus, the model contains only eight material constants excluding the optional relaxation term with p9, otherwise nine constants. Considering the number and complexity of the reproduced effects of the material behavior and the possible fitting quality, verified in the next section, the number of parameters is quite low. Accordingly, parameter optimization is relatively easy, fast and robust.
3. Homogeneous deformations In all following calculations, also the FEM-applications, the model constant p9 is throughout set to zero. Thus, relaxation is not taken into account in any way. Moreover, for all simulations the same set of constants p1 to p8 and K is used. p1=0,035 MPa p2=0,37 MPa p3=0,17 p4=2,4 K =4000 MPa
p5 =0,01 MPa p6=6,4 p7 =5,5 p8=0,24 MPa
Fig. 6. Simulated transient cycles of the tension test, compared with experimental data (dotted lines).
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
Fig. 7. Simulated cyclic tension loading with prescribed load.
Fig. 8. Dependency of small cycles on the prestraining.
1029
1030
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
These values are determined during an optimization process, within which the stationary cycles of the tension- as well as of the shear-test (both Fig. 1) were integrated. Such an optimization strategy is preferred, because the stationary cycles are usually of maximal importance. Fig. 5 shows the comparison of experimental and calculated data of the stationary cycles. Obviously, the agreement is fairly good. Also the transient cycles in Fig. 6 show the correct characteristic, however the agreement to the experimental data is
Fig. 9. Simple shear with rotating axes.
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1031
not as good as that of the stationary cycles. This is not at least a consequence of the fact, that these cycles are not included in the optimization process. A special quality of the constitutive model is the fast transition from the first transient to an almost stationary cycle. However, this is a special effect of the realized kind of loading process. In Fig. 7 another loading process with a constant upper force bound instead
Fig. 10. Cyclic bending of cantilever beam.
1032
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
of an elongation bound is presented. As in experimental data a stationary cycle is reached not until a multitude of previous cycles. In technical applications often periodically changing loadings with relatively small amplitudes around a preloaded mean value are of great interest. It is well-known, that the mean gradient of those small cycles is reduced by an earlier prestraining. Also this effect is reproducible by the constitutive model. Figure 8 shows stationary small cycles with the same elongation bounds after different prestrainings. Obviously, not only the stress values, but also the mean gradient of the cycles are reduced with increasing prestraining. The last example of homogeneous deformations is to be a simple shear-process with rotating axes. This kind of loading is introduced by Gent (1960) and Ahmadi et al. (1999). Rubber exhibits its inelastic behavior by a non-parallelism of shear and force direction in the stationary state. The constitutive model acts in a similar manner. Figure 9 shows the course of the simulated process in the case of a prescribed amount of the force vector and the phase angles between shear and force direction with different reached shear amplitudes. Fortunately, in all these loading situations the constitutive model exhibits a simulated material behavior, which is near to the expected behavior of real rubber.
Fig. 11. Cylinder with spherical hollow under axial pressure.
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1033
4. FEM-applications Beyond the simulations of homogeneous loadings, two inhomogeneous loading examples, calculated with the help of the FEM, shall be presented and discussed. The goal of these intentionally simple examples is to demonstrate the application of the presented constitutive model within the FEM and to discuss some interesting consequences of the typical rubber material behavior rather than to treat current practical questions. The first example is the cyclic bending of a cantilever beam as a plain strain process. The cyclic altering bending force is prescribed. The deformations of the maximum bent states and the following unloaded states of the three cycles in the beginning are displayed in Fig. 10. As expected, the displacements increase in both states from cycle to cycle, while the increases become smaller at each cycle. The remaining deformations after unloading are considerable, in particular in the vicinity of the fixing, where the strongest bending moment appears.
Fig. 12. Stress distribution during maximum deformation. Linear scaling with 13 intervals between 0.78 and 1.07 MPa.
1034
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
The second example is calculated as a rotationally symmetrical problem. A plane of symmetry is modeled with appropriate boundary conditions. The actual simulated structure is a cylinder 40 mm long being 30 mm in diameter with a spherical shaped hollow being 16 mm in diameter in the middle of the cylinder, which is loaded by uniform pressure on the frontal surfaces. Once more the specimen is loaded in a cyclic manner with pressure changes between zero and a constant maximum value. In Fig. 11 the undeformed mesh and the maximum deformed state of the first and second cycle are outlined. As in the above mentioned bending example the deformations increase from cycle to cycle. Figure 12 presents the stresses at the maximum deformation during an almost stationary cycle. In each point the higher one of the two stress eigenvalues belonging to the eigendirections within the output area is plotted. The maximum stress amount corresponds to tension stresses appearing just at the intersection point of the median line with the surface of the hollow. Interestingly, these stresses are slightly increasing during the first cycles, until stationarity is reached (about 6%). Apart from this the
Fig. 13. Remaining stresses after unloading. Linear scaling; 1st cycle: 13 intervals between 0.020 and 0.028 MPa; 2nd cycle: 14 intervals between 0.022 and 0.026 MPa; 3rd cycle: 13 intervals between 0.022 and 0.030 MPa.
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1035
stress distributions of all cycles at maximum load are similar to each other in this example. The stresses in the unloaded state are much smaller than in the deformed state, but they do not vanish because of the inhomogeneous previous deformations. In Fig. 13 the stress distributions of the unloaded states of the first three cycles are shown like the stresses in Fig. 12. The existence of these stress distributions indicates that the material is softened due to the loading history in a varying manner. Thus, to examine the stresses within extensive structures, the consideration of the softening effects in the constitutive model is of vital importance.
5. Conclusions The constitutive model presented reproduces several effects, which are typical for the material behavior of rubber, like strong nonlinearities, hysteresis and softening effects due to the loading history. The final form is completely three-dimensional, all equations contain Lagrangean stress and deformation tensors exclusively. In the incompressible and intrinsic form only eight constants are to be identified. Optional model constants are the bulk modulus for almost incompressible formulations and one other parameter for the integration of a relaxation term. In view of the number and complexity of the reproduced effects these numbers of model constants are quite low. Accordingly, parameter identification is relatively easy. With fitted constants the constitutive model produces a surprisingly good agreement with the stationary cycles of tension and shear tests simultaneously. Simulations of other homogeneous deformations yield expected results, too. Finally, exemplary applications within the Finite Element Method show the suitability and efficiency of the presented constitutive model for the implementation in this numerical method and for practical applications, including realistic predictions of the stress and strain distributions within the material.
References Ahmadi, H.R., Gough, J., Muhr, A.H., 1999. Bi-axial experimental techniques highlighting the limitations of strain-energy description of rubber. In: Dorfmann, A., Muhr, A. (Eds.), Constitutive Models for Rubber. Balkema, Rotterdam, pp. 65–71. Bergstro¨m, J.S., Boyce, M.C., 1998. Constitutive modeling of the large strain time-dependent behavior of elastomers. J. Mech. Phys. Solids 46, 931–954. Gent, A.N., 1960. Simple rotary dynamic testing machine. Br. J. Appl. Phys. 11, 165. Houlsby, G.T., Puzrin, A.M., 2000. A thermomechanical framework for constitutive models for rateindependent dissipative materials. International Journal of Plasticity 16, 1017–1047. Kaliske, M., Rothert, H., 1999. Viscoelastic and elastoplastic damage formulations. In: Dorfmann, A., Muhr, A. (Eds.), Constitutive Models for Rubber. Balkema, Rotterdam, pp. 159–167. Khan, A., Lopez-Pamies, O., 2002. Time and temperature dependent response and relaxation of soft polymer. International Journal of Plasticity, in press. Khan, A., Zhang, H., 2001. Finite deformation of polymer: experiments and modeling. International Journal of Plasticity 17, 1167–1188.
1036
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
Lion, A., 2000. Constitutive modelling in finite thermoviscoplasticity: a physical approach based on nonlinear rheological models. International Journal of Plasticity 16, 469–494. Lion, A., Haupt, P., 2001. Finite viscoelasticity: a physical approach based on relaxation spectra, fractional derivatives and process-dependent viscosities. In: Besdo, D., Schuster, R.H., Ihlemann, J. (Eds.), Constitutive Models for Rubber II. Swets & Zeitlinger, Lisse, pp. 65–72. Lubarda, V.A., Benson, D.J., Meyers, M.A., 2002. Strain-rate effects in rheological models of inelastic response. International Journal of Plasticity, in press. Lulei, F., Miehe, C., 2001. A physically-based constitutive model for finite viscoelastic deformations in rubbery polymers based on directly evaluated micro-macro-transition. In: Besdo, D., Schuster, R.H., Ihlemann, J. (Eds.), Constitutive Models for Rubber II. Swets & Zeitlinger, Lisse, pp. 117–125. Makosey, S.J., Rajagopal, K.R., 2001. The application of ideas associated with materials with memory to modeling the inelastic behavior of solid bodies. International Journal of Plasticity 17, 1087–1117. Puzrin, A.M., Houlsby, G.T., 2001. A thermomechanical framework for rate-independent dissipative materials with internal functions. International Journal of Plasticity 17, 1147–1165. Reese, S., Govindjee, S., 1998. A Theory of finite viscoelasticity and numerical aspects. Int. J. Solids Structures 35, 3455–3482. Sussmann, T., Bathe, K.-J., 1987. A Finite Element Formulation for Nonlinear Incompressible Elastic and Inelastic Analysis. Computers and Structures 26, 357. Truesdell, C., Noll, W., 1965. The non-linear field theories of mechanics. In: Flu¨gge, S. (Ed.), Encyclopedia of Physics. Springer, Berlin, pp. 20–35. Xiao, H., Bruhns, O.T., Meyers, A., 2000. A consistent finite elastoplasticity theory combining additive and multiplicative decomposition of the stretching and the deformation gradient. International Journal of Plasticity 16, 143–177. Xiao, H., Bruhns, O.T., Meyers, A., 2001. Large strain responses of elastic-perfect plasticity and kinematic hardening plasticity with the logarithmic rate: swift effect in torsion. International Journal of Plasticity 17, 211–235.
A phenomenological constitutive model for rubberlike materials and its numerical applications D. Besdo, J. Ihlemann* Institut fu¨r Mechanik, Universita¨t Hannover, Appelstr. 11, 30167 Hannover, Germany
Abstract A new phenomenological inelastic constitutive model for rubberlike materials is presented and its good correspondence to cyclic measurements is demonstrated for both uniaxial tension- and simple shear-tests up to large deformations. The model implies strong nonlinearities, hysteresis, the influence of the loading history and remaining deformations after unloading. Results of finite element calculations are represented to show the suitability of the constitutive model within this method for practical applications. # 2003 Elsevier Science Ltd. All rights reserved. Keywords: Rubber; Constitutive model; Inelasticity; FEM; Large deformation
1. Introduction Numerous technical applications of rubberlike material lead to large deformations. Furthermore, inelastic effects are often of great significance to the function of rubber containing structures. To simulate the mechanical behavior of those continua usually the finite element method is applied, in which by now the purely numerical problems due to rubber properties like strong nonlinearity or vanishing compressibility are controllable fairly well. On the other hand there is still a need for suitable three-dimensional constitutive models. Up to now the molecular foundations of the inelastic behavior of filled vulcanizates under large deformations, so far as known, are not yet transferable to suitable constitutive equations for the FEM in a satisfactory manner.
* Corresponding author. Tel.: +49-511-762-4120; fax: +49-511-762-4182. E-mail address: [email protected] (J. Ihlemann). 0749-6419/03/$ - see front matter # 2003 Elsevier Science Ltd. All rights reserved. PII: S0749-6419(02)00091-8
1020
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
Thus, phenomenological constitutive models are the remaining tools to make simulations possible. To consider the material behavior phenomenologically, suitable and extensive measurements are needed to identify those characteristics which are to be modeled (as an example see Khan and Lopez-Pamies, 2002). Preferably, the boundary conditions of those experiments should be as easy and unique as possible. Homogeneous deformations are preferred to allow the precise consideration of the material behavior itself instead of geometrical influences. The constitutive model, which will be presented in this context, is intended to be useful to the simulation of loading processes with rubber–typical large deformations. Moreover, strong nonlinearities, hysteresis and softening effects caused by the loading history are to be included in the model to allow the consideration of stress distributions within complex structures at the end. For that purpose, special cyclic tension and shear tests with varying upper elongation bounds have turned out to be useful. In Fig. 1 the results of an uniaxial tension test and a symmetrical simple shear test are outlined, both carried out by the Dutch company Akzo Nobel with a carbon black filled rubber used in air springs. Starting each with an initially unloaded specimen, both experiments consist of five increasing upper elongation bounds with several cycles at each level, approaching stationarity. As desired, the interesting characteristics of the material behavior are easy to see.
2. Constitutive model Constitutive descriptions are expressed either in Eulerian or in Lagrangean denotation. Unfortunately, the first Piola–Kirchhoff-stress tensor T, whose coefficients are used in Fig. 1, does not belong to one of these denotations. Hence, a transformation into Cauchy-stresses or second Piola–Kirchhoff-stresses T~ is necessary. For the corresponding transformation of the experimental data, an ideal incompressibility is assumed. The stationary cycles of the tension-test in terms of second Piola–Kirchhoff-stresses are plotted in Fig. 2. Obviously, the curvature of the lines is accentuated. The loading as well as the unloading curve have a changing sign of the curvature. In an Eulerian description (Fig. 3) using Cauchy-stresses no change of the curvature sign can be recognized. Therefore, this description seems to be simpler and, hence, to be preferred for the development of a constitutive model. In the case of the shear-test the curves for the Eulerian description do not differ from the corresponding graph in Fig. 1. However, the final aim will be a formulation using the quantities of the Lagrangean denotation, since most finite element codes use a Lagrangean formulation of the constitutive equations for implementation. Moreover, in order to implement the constitutive model, of course it must be formulated in a general tensorial 3D-scheme. A general problem of simulating rubber behavior is the near-incompressibility of those materials. Typically, the ratio of the bulk modulus to the shear modulus is several thousand. Because of this, the material behavior concerning pure changes in volume and on the other hand concerning deformations without any density
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1021
Fig. 1. Results of cyclic uniaxial tension and symmetric simple shear tests carried out by the Dutch company Akzo Nobel (transient cycles are plotted only partially). First Piola–Kirchhoff-stress tensor T vs. deformation gradient tensor F.
1022
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
changes preferably should be modeled and calculated separately. The volume ratio and thus, the ratio of the reference density to the current density, is connected with J3, the third main invariant of the deformation gradient F, which is equal to the determinant of the matrix of the coefficients of F in any orthonormalized coordinates. J3 ¼ I3 ðFÞ ¼ det½Fab With the help of J3 modifications of usual quantities are calculable, which are insensitive to the magnitude of the volume ratio (Sussmann and Bathe, 1987). Those modified quantities will be characterized by a superscript M. The use of these modified quantities corresponds to the idea of the ‘surface of natural density’ described by Makosey and Rajagopal (2001). In this context, corresponding modifications of
the right Cauchy–Green tensor C and the Jaumann rate b of the left Cauchy–Green tensor b are needed (0 denotes the purely mathematical deviatoric part). M
C¼
J32=3 C;
0 M 1 b ¼ b b b;
M
M
b ¼ J32=3 b
The new constitutive model should be able to describe the stationary cycles and the transition to higher strain levels as precisely as possible without a too high number of material constants. Nevertheless, the hysteresis, the remaining strains
Fig. 2. Stationary cycles of the tension test in Lagrangean description. Second Piola–Kirchhoff-stress tensor T~ vs. Green–Lagrange strain tensor .
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1023
after unloading and the influence of the strain history should be modeled in satisfactory approximation. The strain history is fitted by a scalar parameter. Thus, the influence of old loading directions leading to a material anisotropy is completely neglected. To define the one parameter, a special invariant is used, which represents the maximum difference of two eigenvalues. XT ¼ maxfjXI XII j; jXII XIII j; jXIII XI jg Applied to the Cauchy-stresses this is the reduced stress corresponding to Tresca’s yield criterion. The scalar parameter representing the strain history, is to be the maximum value of this invariant of the volume ratio insensitive modifications of the Cauchy–Green tensors in the past. M M H CT ðtÞ ¼ max CT ð Þ; 0 4 4 t M
In the following CH T is called the history function. Several constitutive models for rubber have been developed using rheological models with serial and parallel combinations of elastic, viscous and plastic elements
Fig. 3. Stationary cycles of the tension test in Eulerian description. Cauchy stress tensor vs. left Cauchy–Green tensor b.
1024
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
(Bergstro¨m and Boyce, 1998; Kaliske and Rothert, 1999; Lion, 2000; Lion and Haupt, 2001; Lubarda et al., 2002; Lulei and Miehe, 2001; Reese and Govindjee, 1998). Instead of this, the formulation of the present paper is in general based on the fundamental idea, to approximate the measured cycles by two limiting lines as asymptotic curves. The changes of the stresses in the range prescribed by these limits, which will actually be limiting tensors in the model, may be described by a specific differential equation, while strictly speaking the complete stresses will be an additive composition of some basic stresses and the solution of this differential equation (see Fig. 4). Additional experiments with varying strain rates showed that although the succession of the deformations has a strong influence on the material behavior, the influence of the amount of the velocity is very small over a wide velocity range. Therefore, the material description is designed to be independent of the rate in the above mentioned sense. Due to this feature this constitutive model is called an intrinsic model. Those rate independent models are often used to describe the behavior of metals undergoing finite deformations (Houlsby and Puzrin, 2000; Makosey and Rajagopal, 2001; Puzrin and Houlsby, 2001), whereas many other materials, e.g. teflon (Khan and Zhang, 2001), exhibit a distinctly time dependent behavior.
Fig. 4. Computed stationary tension cycle with corresponding basic and limiting stresses.
1025
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
The stresses coming from the solution of the differential equation, providing the determinant part of the hysteresis, are identified by means of the (additional) second Piola–Kirchhoff-stresses T~ A . The remaining parts of the total stresses are the so called basic stresses, which are similar to neo Hookean stresses but not identical with them, because the parameter is not constant in this case, but it is a definite function of three material constants p1, p2, p3 and in particular the history function M
CH T ðtÞ. Altogether, the stresses are divided into one part, which is physically corresponding to the deviatoric part of the Cauchy stress tensor and another part, representing the hydrostatic material behavior containing the bulk modulus K. 0 M A ~ ~ T ¼ 2C þ T C C1 þ K J3 ðJ3 1ÞC1 with :
p2 ¼ p1 þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi M
2
1 þ p3 CH T
Supposing a totally incompressible material, the bulk modulus tends to infinity and the value of the volume ratio J3 is identical with one in any case. Thus, the hydrostatic pressure is indeterminate and has to be calculated directly from the equilibrium equations. The differential equation for the additional stresses T~ A should cause these stresses A ~ T to tend to the current limiting stresses T~ L . This requirement leads to an
expression for T~ A , the simple Lagrangean rate of T~ A . T~ A ¼ LT T~ L T~ A
M
is again a definite function of material constants and the history function CH T ðtÞ (see below). LT is an invariant of a Lagrangean equivalent to an Eulerian time derivative preserving the required intrinsic character of the model. The tensor L will
be introduced below. The combination of the time derivative T~ A and the time dependent invariant LT cause a similar effect as the use of the ‘inelastic strain pathlength’ of Makosey and Rajagopal (2001). Beyond the frame of the intrinsic description a relaxation term can be easily inserted if required. T~ A ¼ LT T~ L T~ A p9 T~ A
For the final version of the differential equation, an additional term is needed to prevent physically hydrostatic parts of the limiting stresses influencing the evolution of the physically deviatoric parts of the stresses T~ A .
1026
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036 1 T~ A ¼ LT T~ L T~ A p9 T~ A þ T~ A C C1 3 p4 with : ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi M 1 þ p3 CH T
As mentioned above, the formulation of the limiting stresses will be easier in the Eulerian description. These stresses are explained as a tensorial function (Truesdell and Noll, 1965) of the Jaumann rate of the left Cauchy–Green tensor (for a discussion of useful objective corotational rates see Xiao et al., 2000, 2001). To retain the intrinsic character of the description, this time derivative is divided by its own M
invariant b T . Moreover, the vicinity to the prestraining maximum is included with M
the help of the ratio of the current value of b T with its maximum value in the past, which is equal to the history function. 2 L ¼
00
M M
11
M
3
BB CC 7 16 6 expBBp7 b b T CC þ p8 b 7 @ @ A A 4 5 M M M J3 bH T bT bT
M
with :
M
b T ¼ CT
To transfer Mthis equation identically into the Lagrangean notation, a Lagrangean M
equivalent to b is introduced, which has the same eigenvalues as b itself. However, it should be mentioned, that this tensor L may be unsymmetric, but it is still suitable to be an argument of a tensorial function. 0 M 1 1 1 C C þ C C L¼ C 2 With the help of L the transfer to a pure Lagrangean notation is possible without further problems. Once more, is a function of material constants and the history M
function CH T ðtÞ.
2
3 11 M L C L TA A þ p8 5 C1 T~ L ¼ 4 exp@@p7 LT M H LT CT 0 with :
00
1
B C M B C p6 B1 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C ¼ p5 C H TB 2 C @ A M p26 þ CH T
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
Fig. 5. Tension and shear test: comparison of experimental and calculated data.
1027
1028
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
With the explanation of the limiting stresses the constitutive model is completed and ready for material behavior simulation and implementation in the finite element method. Apart from the bulk modulus, the model contains only eight material constants excluding the optional relaxation term with p9, otherwise nine constants. Considering the number and complexity of the reproduced effects of the material behavior and the possible fitting quality, verified in the next section, the number of parameters is quite low. Accordingly, parameter optimization is relatively easy, fast and robust.
3. Homogeneous deformations In all following calculations, also the FEM-applications, the model constant p9 is throughout set to zero. Thus, relaxation is not taken into account in any way. Moreover, for all simulations the same set of constants p1 to p8 and K is used. p1=0,035 MPa p2=0,37 MPa p3=0,17 p4=2,4 K =4000 MPa
p5 =0,01 MPa p6=6,4 p7 =5,5 p8=0,24 MPa
Fig. 6. Simulated transient cycles of the tension test, compared with experimental data (dotted lines).
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
Fig. 7. Simulated cyclic tension loading with prescribed load.
Fig. 8. Dependency of small cycles on the prestraining.
1029
1030
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
These values are determined during an optimization process, within which the stationary cycles of the tension- as well as of the shear-test (both Fig. 1) were integrated. Such an optimization strategy is preferred, because the stationary cycles are usually of maximal importance. Fig. 5 shows the comparison of experimental and calculated data of the stationary cycles. Obviously, the agreement is fairly good. Also the transient cycles in Fig. 6 show the correct characteristic, however the agreement to the experimental data is
Fig. 9. Simple shear with rotating axes.
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1031
not as good as that of the stationary cycles. This is not at least a consequence of the fact, that these cycles are not included in the optimization process. A special quality of the constitutive model is the fast transition from the first transient to an almost stationary cycle. However, this is a special effect of the realized kind of loading process. In Fig. 7 another loading process with a constant upper force bound instead
Fig. 10. Cyclic bending of cantilever beam.
1032
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
of an elongation bound is presented. As in experimental data a stationary cycle is reached not until a multitude of previous cycles. In technical applications often periodically changing loadings with relatively small amplitudes around a preloaded mean value are of great interest. It is well-known, that the mean gradient of those small cycles is reduced by an earlier prestraining. Also this effect is reproducible by the constitutive model. Figure 8 shows stationary small cycles with the same elongation bounds after different prestrainings. Obviously, not only the stress values, but also the mean gradient of the cycles are reduced with increasing prestraining. The last example of homogeneous deformations is to be a simple shear-process with rotating axes. This kind of loading is introduced by Gent (1960) and Ahmadi et al. (1999). Rubber exhibits its inelastic behavior by a non-parallelism of shear and force direction in the stationary state. The constitutive model acts in a similar manner. Figure 9 shows the course of the simulated process in the case of a prescribed amount of the force vector and the phase angles between shear and force direction with different reached shear amplitudes. Fortunately, in all these loading situations the constitutive model exhibits a simulated material behavior, which is near to the expected behavior of real rubber.
Fig. 11. Cylinder with spherical hollow under axial pressure.
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1033
4. FEM-applications Beyond the simulations of homogeneous loadings, two inhomogeneous loading examples, calculated with the help of the FEM, shall be presented and discussed. The goal of these intentionally simple examples is to demonstrate the application of the presented constitutive model within the FEM and to discuss some interesting consequences of the typical rubber material behavior rather than to treat current practical questions. The first example is the cyclic bending of a cantilever beam as a plain strain process. The cyclic altering bending force is prescribed. The deformations of the maximum bent states and the following unloaded states of the three cycles in the beginning are displayed in Fig. 10. As expected, the displacements increase in both states from cycle to cycle, while the increases become smaller at each cycle. The remaining deformations after unloading are considerable, in particular in the vicinity of the fixing, where the strongest bending moment appears.
Fig. 12. Stress distribution during maximum deformation. Linear scaling with 13 intervals between 0.78 and 1.07 MPa.
1034
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
The second example is calculated as a rotationally symmetrical problem. A plane of symmetry is modeled with appropriate boundary conditions. The actual simulated structure is a cylinder 40 mm long being 30 mm in diameter with a spherical shaped hollow being 16 mm in diameter in the middle of the cylinder, which is loaded by uniform pressure on the frontal surfaces. Once more the specimen is loaded in a cyclic manner with pressure changes between zero and a constant maximum value. In Fig. 11 the undeformed mesh and the maximum deformed state of the first and second cycle are outlined. As in the above mentioned bending example the deformations increase from cycle to cycle. Figure 12 presents the stresses at the maximum deformation during an almost stationary cycle. In each point the higher one of the two stress eigenvalues belonging to the eigendirections within the output area is plotted. The maximum stress amount corresponds to tension stresses appearing just at the intersection point of the median line with the surface of the hollow. Interestingly, these stresses are slightly increasing during the first cycles, until stationarity is reached (about 6%). Apart from this the
Fig. 13. Remaining stresses after unloading. Linear scaling; 1st cycle: 13 intervals between 0.020 and 0.028 MPa; 2nd cycle: 14 intervals between 0.022 and 0.026 MPa; 3rd cycle: 13 intervals between 0.022 and 0.030 MPa.
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
1035
stress distributions of all cycles at maximum load are similar to each other in this example. The stresses in the unloaded state are much smaller than in the deformed state, but they do not vanish because of the inhomogeneous previous deformations. In Fig. 13 the stress distributions of the unloaded states of the first three cycles are shown like the stresses in Fig. 12. The existence of these stress distributions indicates that the material is softened due to the loading history in a varying manner. Thus, to examine the stresses within extensive structures, the consideration of the softening effects in the constitutive model is of vital importance.
5. Conclusions The constitutive model presented reproduces several effects, which are typical for the material behavior of rubber, like strong nonlinearities, hysteresis and softening effects due to the loading history. The final form is completely three-dimensional, all equations contain Lagrangean stress and deformation tensors exclusively. In the incompressible and intrinsic form only eight constants are to be identified. Optional model constants are the bulk modulus for almost incompressible formulations and one other parameter for the integration of a relaxation term. In view of the number and complexity of the reproduced effects these numbers of model constants are quite low. Accordingly, parameter identification is relatively easy. With fitted constants the constitutive model produces a surprisingly good agreement with the stationary cycles of tension and shear tests simultaneously. Simulations of other homogeneous deformations yield expected results, too. Finally, exemplary applications within the Finite Element Method show the suitability and efficiency of the presented constitutive model for the implementation in this numerical method and for practical applications, including realistic predictions of the stress and strain distributions within the material.
References Ahmadi, H.R., Gough, J., Muhr, A.H., 1999. Bi-axial experimental techniques highlighting the limitations of strain-energy description of rubber. In: Dorfmann, A., Muhr, A. (Eds.), Constitutive Models for Rubber. Balkema, Rotterdam, pp. 65–71. Bergstro¨m, J.S., Boyce, M.C., 1998. Constitutive modeling of the large strain time-dependent behavior of elastomers. J. Mech. Phys. Solids 46, 931–954. Gent, A.N., 1960. Simple rotary dynamic testing machine. Br. J. Appl. Phys. 11, 165. Houlsby, G.T., Puzrin, A.M., 2000. A thermomechanical framework for constitutive models for rateindependent dissipative materials. International Journal of Plasticity 16, 1017–1047. Kaliske, M., Rothert, H., 1999. Viscoelastic and elastoplastic damage formulations. In: Dorfmann, A., Muhr, A. (Eds.), Constitutive Models for Rubber. Balkema, Rotterdam, pp. 159–167. Khan, A., Lopez-Pamies, O., 2002. Time and temperature dependent response and relaxation of soft polymer. International Journal of Plasticity, in press. Khan, A., Zhang, H., 2001. Finite deformation of polymer: experiments and modeling. International Journal of Plasticity 17, 1167–1188.
1036
D. Besdo, J. Ihlemann / International Journal of Plasticity 19 (2003) 1019–1036
Lion, A., 2000. Constitutive modelling in finite thermoviscoplasticity: a physical approach based on nonlinear rheological models. International Journal of Plasticity 16, 469–494. Lion, A., Haupt, P., 2001. Finite viscoelasticity: a physical approach based on relaxation spectra, fractional derivatives and process-dependent viscosities. In: Besdo, D., Schuster, R.H., Ihlemann, J. (Eds.), Constitutive Models for Rubber II. Swets & Zeitlinger, Lisse, pp. 65–72. Lubarda, V.A., Benson, D.J., Meyers, M.A., 2002. Strain-rate effects in rheological models of inelastic response. International Journal of Plasticity, in press. Lulei, F., Miehe, C., 2001. A physically-based constitutive model for finite viscoelastic deformations in rubbery polymers based on directly evaluated micro-macro-transition. In: Besdo, D., Schuster, R.H., Ihlemann, J. (Eds.), Constitutive Models for Rubber II. Swets & Zeitlinger, Lisse, pp. 117–125. Makosey, S.J., Rajagopal, K.R., 2001. The application of ideas associated with materials with memory to modeling the inelastic behavior of solid bodies. International Journal of Plasticity 17, 1087–1117. Puzrin, A.M., Houlsby, G.T., 2001. A thermomechanical framework for rate-independent dissipative materials with internal functions. International Journal of Plasticity 17, 1147–1165. Reese, S., Govindjee, S., 1998. A Theory of finite viscoelasticity and numerical aspects. Int. J. Solids Structures 35, 3455–3482. Sussmann, T., Bathe, K.-J., 1987. A Finite Element Formulation for Nonlinear Incompressible Elastic and Inelastic Analysis. Computers and Structures 26, 357. Truesdell, C., Noll, W., 1965. The non-linear field theories of mechanics. In: Flu¨gge, S. (Ed.), Encyclopedia of Physics. Springer, Berlin, pp. 20–35. Xiao, H., Bruhns, O.T., Meyers, A., 2000. A consistent finite elastoplasticity theory combining additive and multiplicative decomposition of the stretching and the deformation gradient. International Journal of Plasticity 16, 143–177. Xiao, H., Bruhns, O.T., Meyers, A., 2001. Large strain responses of elastic-perfect plasticity and kinematic hardening plasticity with the logarithmic rate: swift effect in torsion. International Journal of Plasticity 17, 211–235.