Computational modelling of rubber-like materials under monotonic and cyclic loading

(1)

(2)

(3)

Computational Modelling of Rubber-like Materials

under Monotonic and Cyclic Loading

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof.dr.ir. J.T. Fokkema, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op vrijdag 10 november 2006 om 10:00 uur

door

Zhanqi GUO

(4)

Dit proefschrift is goedgekeurd door de promotor:

Prof.dr.ir. L.J. Sluys

Samenstelling promotiecommissie

Rector Magnificus Voorzitter

Prof.dr.ir. L.J. Sluys Technische Universiteit Delft, promotor Prof. R.W. Ogden University of Glasgow, UK

Prof. dr.ir. E. van der Giessen Rijksuniversiteit Groningen Prof.dr.ir. K.van Breugel Technische Universiteit Delft Prof.dr.ir. P. Stroeven Technische Universiteit Delft

Hoogleraar: Northern Jiaotong University, China Dr. A. Scarpas Technische Universiteit Delft

Dr.ir. G.N. Wells Technische Universiteit Delft

Prof.dr.ir. J. Blaauwendraad Technische Universiteit Delft, reservelid

Computational Modelling of Rubber-like Materials under Monotonic and Cyclic Loading / Z. Guo

Thesis Delft University of Technology. With ref. – With summary in Dutch. ISBN-10: 90-9021154-3

ISBN-13: 978-90-9021154-1

Keywords: Rubber-like material, Constitutive modelling, Mullins effect, Continuum damage mechanics, Pseudo-elastic model, Finite element method Cover design: Jian Guo

Copyright @2006 by Zhanqi Guo

All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording or by any information storage and retrieval system, without written permission from the copyright owner.

(5)

v

Acknowledgements

The research presented in this dissertation has been carried out in the framework of a doctoral study at Delft University of Technology, Faculty of Civil Engineering and Geosciences.

I would like to express my sincere gratitude to my supervisor Prof.dr.ir. L.J. Sluys for his scientific and mental encouragement, continuous support and valuable guidance in the course of my study. Without his encouragement and generous help, I would not have been able to complete this thesis.

I am very grateful to Prof.dr.ir. P. Stroeven. I am deeply moved by his active involvement in international science cooperation, his enthusiastic approach to scientific research and optimistic attitude towards life. I cherish our friendship, and highly appreciate the interesting social communications and his patient help in all aspects.

I would like to thank the late Prof.dr.ir. Ch.F. Hendriks, who gave me the opportunity to study and work at Delft University of Technology. I cannot forget his support and help. I would also like to show my thanks to Prof.dr.ir. K.van Breugel for his support and help after he took over the position of head of the material science group.

Moreover, I would like to acknowledge all staff members of the Section of Materials Science and Sustainable Construction, and the Section of Structural Mechanics, as well as my Chinese friends in Delft for their cooperation and support. Sincere gratitude is also given to the CICAT staff members, Drs. P. Althuis and C. Timmers for support and help.

My special thanks go to Dr. Gregory Chagnon and his colleagues at the French Research Department for providing the experimental data, B.B.C. Jong and M.A.N. Hendriks from TNO DIANA for upgrading the DIANA user-subroutine, and Frank Everdij and Frank Custers for maintaining my computer system.

I cannot end without thanking my family for their long-distance encouragement and support throughout my overseas time. Last but not least, I am very thankful to my wife, Xiuyun, and my son, Jian, for their lasting love, support and patience.

Zhanqi Guo

(6)

(7)

vii

Chapter 1 Introduction

1.1 Background

The word “rubber” originally comes from the natural product of the tree Hevea Braziliensis, of which the chief chemical constituent is the rubber hydrocarbon (C5H8)n or

polyisoprene. The only other rubber then known was gutta-percha or balata, another natural product represented by the same empirical formula (C5H8)n, but differing slightly in the

structural form of the molecule. Recent years have seen the development of a very large number of synthetic rubbers having a wide variety of chemical constitutions. To further improve the stiffness and the strength of the rubber a variety of additive fillers is also employed. The term “rubber” in my studies is now employed to include both natural rubber and synthetic rubbers. Natural rubber as well as synthetic rubber falls in the class of rubber-like materials.

Rubber-like materials exhibit a highly nonlinear behaviour characterized by hyper-elastic deformability and incompressibility or near-incompressibility. Normally, the maximum

extensibility of rubber could reach values varying from 500% to 1000% and the typical stress-strain curve in tension is markedly nonlinear so that Hooke’s law cannot be applied and it is not possible to assign a definite value to the Young’s modulus except in the region of small strains, where the Young’s modulus is of the order of 1MPa. In contrast, the Young’s modulus for typical hard solids is in the region 104_-106_{MPa and the maximum elastic}

extensibility of hard solids seldom exceeds 1%. Rubber-like materials are effectively incompressible in most situations. However, all real materials are compressible to a certain degree even if the bulk modulus of the rubber-like materials is several orders of magnitude larger than the shear modulus.

In addition, when a rubber specimen is subjected to cyclic loading, the Mullins effect or

stress softening, which is characterized by an important loss of stiffness during the first few cycles, has been observed. Moreover, a carbon black filled rubber after loading and subsequent unloading, in general, does not return to its initial state corresponding to the natural stress-free configuration. But, a residual strain or permanent deformation remains. Finally, the mechanical behaviour of rubber-like materials is time-dependent which can be demonstrated by relaxation and creep experiments.

(12)

2 Chapter 1

An increase of applications requires a better understanding of the mechanical behaviour of rubber-like materials. There is an enormous difference between rubbers and ordinary hard solids. Unlike metals, which require relatively few properties to characterize their behaviour, the behaviour of rubber is complex. Usually, the mechanical behaviour of rubber-like materials can not be described by a simple stress-strain relation, but by the strain energy function.

Therefore, an important problem in non-linear elasticity theory is to come up with a reasonable and applicable elastic law (strain energy function), which is the key to the development of reliable analysis tools. Many attempts have been made to develop a theoretical stress-strain relation that can fit experimental results for hyperelastic materials. For example, Mooney (1940), Rivlin (1948b), Treloar (1944, 1975), Knowles & Sternberg (1973), Ogden (1972a, b) and others have proposed strain energy functions for rubber-like materials (more details can be found in Chapter 4 and relevant references). An excellent agreement has been obtained between Ogden’s formula and Treloar’s experimental data for extensions of unfilled natural rubber up to 700% (Ogden, 1972a, b). But, it is not convenient because many parameters must be determined when we use the Ogden model. However, the models with a lower number of material parameters often fail in the task of describing the response of a rubber material under different states of deformation without changing the model parameters.

Another significant phenomenon of mechanical behaviour for rubber-like materials is stress softening. Many researches are studying this so-called Mullins effect by means of either phenomenological-based models or molecular-based models. Phenomenological models are used to describe the Mullins effect based purely on experimental observations and molecular-based models are used to explain the effect by the macroscopical nature of the rubber molecular structure. Aiming at engineering application by minimizing the model parameters with proper prediction of the mechanical behaviour requires further study for this phenomenon.

In general a carbon black reinforced rubber after loading and subsequent unloading does not return to its initial state corresponding to the natural stress-free configuration, but exhibits a residual strain or permanent deformation. The permanent deformation combined with stress-softening effects in rubber-like materials results in complex mechanical behaviour and modelling is still at an early stage.

Fortunately, the rapid development of computers and computational tools, such as the finite element method have provided the powerful means to further study, understand and optimize rubber in engineering applications.

1.2 Objectives and scope of this study

The present work was motivated by the need for an appropriate material model, to be used for hyperelastic rubber-like materials, which contains a small number of material parameters and is able to describe the material response for different deformation modes under reasonable high deformation levels.

(13)

Chapter 1 3

phenomenological approach to rubber elasticity, which is not based on molecular or structural concepts but on purely mathematical reasoning and does not explicitly include the physical connection with the underlying mechanisms of deformation. An elastic model proposed by Gao (1997) contains only two model parameters and is used for describing a general, three-dimensional state of deformation under static loading conditions. This constitutive relation was successfully used for theoretical analyses (Gao 2001, 1999, 1997), without experimental verification. Furthermore, an analysis of more complicated deformation with arbitrary geometry cannot be done theoretically and must be analysed by means of numerical methods. For the numerical analysis the finite element method was chosen. The capabilities of this computational model are firstly investigated by basic numerical calculations and discussions on model parameters and Drucker’s stability postulate. Further evaluations of this model have been done by a comparison between numerical prediction and experimental data under different deformations as well as a comparison with other commonly used material models in this study.

The second purpose of this study is concerned with constitutive models, which incorporate stress-softening phenomena. The continuum damage mechanics concept and a pseudo-elastic model combined with Gao’s elastic law are proposed to describe ideal Mullins effect, respectively. These combinations aim at engineering application in terms of minimizing the number of parameters and reasonably describing the Mullins effect. Numerical simulations for different deformations under cyclic loading compared with experimental data and analyses of inhomogeneous problems are presented to demonstrate the capabilities of the two models.

Finally, the mechanical behaviour of rubber-like materials contains stress softening with residual deformations. A specific constitutive model to capture the Mullins effect and its corresponding permanent deformation is proposed. Comparisons between the numerical simulations and experimental results are used for calibration of the proposed model.

1.3 Outline of the thesis

Chapter 2 starts with the essential aspects of continuum mechanics, the definitions of the strain measures and basic kinematic relations for the description of large strain analyses.

In Chapter 3 the finite element method applied to rubber-like materials, in which displacements and the pressure-field are separately interpolated, is introduced. An interface framework for implementation of new material models is also illustrated in this part.

Chapter 4 presents an overview of phenomenological models for rubber-like materials, the proposed model and the implementation of the formulation into a finite element code. A discussion on model parameters and Drucker’s stability postulate as well as 2D and 3D element tests under different loading modes is included in order to give insight into this model.

(14)

4 Chapter 1

In Chapter 6 the ideal Mullins effect or stress softening is firstly described. The continuum damage mechanics concept and the pseudo elastic model combined with Gao’s elastic law are used to simulate the ideal Mullins effect. Necessary formulations are derived for implementation of the models into a finite element code. Theoretical and numerical analysis of pure shear deformation illustrates the fundamental characteristics of the proposed models.

Verification and applications of the continuum damage mechanics model and the pseudo elastic model are presented in Chapter 7. Two sets of experimental data under different loads are employed to evaluate the two models. Finally, two applications of inhomogeneous analyses demonstrate the characteristics of the continuum damage mechanics model and the pseudo elastic model. Discussion on the model parameters of the pseudo elastic model gives further understanding of the models.

In Chapter 8, a specific model for stress softening with permanent deformation is proposed. Loading, unloading, reloading and secondary unloading paths are described by different expressions of the damage parameter, which totally involves five model parameters and could be estimated separately according to different branches of evolution curves. Finally, two experiments verify this specific model.

(15)

Chapter 2 Basic notations of hyper-elasticity

Any sample of a solid consists of “particles”. These particles may be in the form of atoms, molecules, long chains of molecules in polymers and crystallites in metals if we observe them at microscopic scale. At macroscopic scale, these particles can not be distinguished and are replaced by a continuous medium. Therefore, a mathematical description of the behaviour of an actual solid can be obtained in terms of continuum parameters. This branch of mechanics is called continuum mechanics. In this chapter, a short recapitulation of continuum mechanics is firstly presented. The basic notation, which is used throughout this thesis and the fundamental description of hyper-elastic materials is given in this chapter.

2.1 Finite strain kinematics

In this section, basic principles of finite deformation kinematics are summarized. The objective is to introduce some fundamental kinematic measures and some basic laws used in the following sections. More comprehensive continuum kinematics overviews can be found in many kinematics textbooks, e.g. continuum mechanics (Chandrasekharaiah and Debnath, 1994).

2.1.1 Deformation gradient tensor

Let us consider an arbitrarily selected configuration C0. We refer to C0 as an initial, reference or also as undeformed configuration. In terms of C0 the body is defined to consist of a set of material points, which lie in one to one correspondence with the set of position vectors x0_{. By definition the correspondence between particle and}_x0_{is bijective so that each}

particle is uniquely labelled by its coordinate vector x0_.

Consider now the actual configuration Ct occupied by the body at time t. The particle labelled by x0 moves to a position x, which depends on x0 and t. One possible description of the motion of the body is given by the vector equation

( )

0_,t

=

x x x (2.1) which defines the material or reference description of the motion of a continuum. Conventionally, the material description is referred to as the Langrangian description.

(16)

6 Chapter 2 a x0 x P0 P u dx0 dx C0 Ct

Fig. 2.1 Material and spatial coordinates.

( )

0 ₌ 0 _,t

x x x (2.2)

exists. When Eq. (2.2) is used as the material description, we name it Eulerian description. The position vectors x0_and_{x are defined as follows:}

0 _;

i i xi i

ξ

= =

x i x i (2.3)

let us consider the infinitesimal material vectors dx and dx0_{which can be related by a second}

order tensor F,

0

dx F x = d (2.4)

The tensor F is called the deformation gradient tensor and can be written as

1 1 1 1 2 3 2 2 2 0 1 2 3 3 3 3 1 2 3 x x x x x x x x x ξ ξ ξ ξ ξ ξ ξ ξ ξ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ x F x (2.5) or in component form: , k k i i x F ξ ∂ = ∂ (2.6)

We refer to displacement vector of the particle as u, we have

0

= −

(17)

Chapter 2 7

then, it is easy to obtain

i ij ij j u F δ ξ ∂ = + ∂ (2.8)

where δij is the Kronecker delta. Expression (2.8) can be in tensor form: ∂

= + ∂

u

F I (2.9)

where I is identity tensor.

2.1.2 Volume changes, Jacobian matrix

Let us consider an infinitesimal body. The volume of the cube is represented by dV0 in the

undeformed configuration and dV in the deformed configuration, respectively. We have an undeformed volume

(

)

0 1 1 2 2 3 3 1 2 3

dV = ∂ξi × ∂ξ i ∂ξ i = ∂ ∂ ∂ξ ξ ξ (2.10) and a deformed volume

(

1 1 2 2

)

3 3

(

1 2 3

)

1 2 3

dV = F∂ × ∂ξ F ξ F∂ =ξ F F F× ∂ ∂ ∂ξ ξ ξ (2.11) Substituting the tangent vector into this equation yields

0

dV =JdV (2.12)

where J is called the Jacobian of the transformation. J represents the ratio of the actual volume of the infinitesimal body to the initial volume and can be calculated by the determinant of the deformation gradient.

( )

1 1 1 1 2 3 2 2 2 1 2 3 3 3 3 1 2 3 det x x x x x x J x x x ξ ξ ξ ξ ξ ξ ξ ξ ξ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ F (2.13)

From the mass conservation law

0 0 dVρ=dVρ (2.14) It yields 0 0 dV J dV ρ ρ = = (2.15)

where ρ0 is the material density of volume dV0 in undeformed configuration and ρ is the

(18)

8 Chapter 2

2.1.3 Strain tensor

Based on the physical argument that the volume of any element with non-zero initial volume cannot have a zero or negative volume after deformation, it holds that

( )

det F >0 (2.16)

Thus, the deformation gradient F has its unique inverse F-1_{. The deformation gradient}_{F gives}

a complete description of the motion of a material particle. We can decompose F into a pure

rotation R and a pure deformation U or V in a multiplicative fashion:

=

F RU (2.17)

=

F VR (2.18)

where the orthogonal tensor R is called rotation tensor and the positive definite tensors U and V are called right and left stretch tensors, respectively. This decomposition is commonly

named the polar decomposition of the deformation gradient (Fig. 2.2).

A strain tensor can be related to the change of length of a line element in the deformed state relative to the reference state. The Green-Lagrange strain tensor can be defined based on the difference of the square of the line element.

(

)

(

)

(

)

2 2 1 1 2 2 2 2 T T T T T T T d d d d d d d d d d d d d d − = − = − = − = − X X X F F F F I C I (2.19)

where the tensor C is referred to as the right Cauchy-Green strain tensor, which is defined as

(19)

Chapter 2 9

T

=

C F F (2.20)

and the tensor

T

=

B FF (2.21)

is called the Finger deformation tensor. Tensor C is related to the undeformed configuration

and tensor B is related to the deformed configuration. Furthermore, the Green-Lagrange

strain tensor is defined as

(

)

(

)

1 1

2 2

T

= F F I− = C I − (2.22)

In terms of the derivatives of the displacement field, substituting Eq. (2.9) into Eq. (2.22) shows that the Green-Lagrange strain tensor can be written as

1 1 2 2 T T ∂ ∂ ∂ ∂ = + + ∂ ∂ ∂ ∂ u u u u _(2.23) or in index notation: 1 1 2 2 T T j i k k ij j i i j u u u u γ ξ ξ ξ ξ ∂ ∂ ∂ ∂ = + + ∂ ∂ ∂ ∂ (2.24)

in the first term of the right hand side the linear kinematic relation is recognized, while the second term is the extension for geometrical nonlinearity.

The right Cauchy-Green strain tensor C is a more useful strain measure than the

Green-Lagrange strain tensor γγγγ for the description of constitutive models of hyperelastic materials. The tensor γγγγ and tensor C are symmetric tensors in contrast to the deformation gradient tensor

F. The right Cauchy-Green strain tensor C contains some tensor invariants that are often used

to formulate the constitutive behaviour, such as

(20)

10 Chapter 2 * 1 2 1 * 3 : I I I − − =C I= (2.30)

In which the operator : is the double dot product, according to : =A Bij ji

A B (2.31) The first and second order derivatives of the invariants with respect to the right-Green strain tensor C will be often used in the expression for the total and incremental stress-strain

relation. From the above expressions we can easily derive that

(21)

Chapter 2 11 33 22 23 33 11 31 22 11 12 2 * 3 ₃₃ ₃₁ ₂₃ 2 12 31 11 12 23 12 23 12 22 31 0 0 0 0 0 0 0 0 0 0 0 ₂ ₂ ₂ 0 ₂ ₂ ₂ 0 0 ₂ ₂ ₂ C C C C C C C C C I _C _C _C C C C C C C C C C C − − − ∂ ₌ ₋ − ∂ − − − − − C (2.40)

Another measure of deformation is the stretch ratio λ, which is the ratio of the length of the material line after and before deformation,

T T T T T T T d d d d d d d d _{d d} _{d d} _{d d} λ= x = x x = F F == C (2.41)

The stretch ratios in the principal directions are principal stretches, which are denoted as

λ1, λ2 and λ3. Thus, the relative volume change can be written as

1 2 3

J =λ λ λ (2.42)

2.2 Stress tensor

The concept of stress is defined as the force per unit load-carrying area. When the displacement gradients remain small compared to unity, it is not relevant in which configuration the load-carrying area is measured, the deformed or the undeformed configuration. However, the matter becomes important if this assumption is no longer applicable. Then, it must be defined unambiguously to which configuration the stresses are referred to. For engineering purposes the Cauchy stress tensor σσσσ is defined to describe the magnitude of the stresses in the current configuration. The Cauchy stress tensor that contains “ true stress components” is equal to the force per unit of deformed area and is the essential stress tensor of continuum mechanics. It satisfies the usual equilibrium equations in the current configuration Ct

div +ρg=0 in C_t (2.43) div is the divergence operator, g is the vector of the acceleration of gravity force. The Cauchy

stress tensor is symmetric

T

= (2.44)

(22)

12 Chapter 2

1 T

J

= F F (2.45)

J is the relative volume change.

2.3 Constitutive relation

The second law of thermodynamics requires that the internal production of entropy must be positive, so that heating can never be less than the rate of heat supply. The expression of this law is referred to as the Clausius-Duhem inequality. For the case of purely mechanical theory, the Clausius-Duhem inequality has the form

: 0

ψ

− + Σ E≥ (2.46)

where Σ is stress, which is conjugate to the particular strain tensor E and ψ denotes the free energy. If ψ is only a function of deformation F, it is also referred to as the strain energy

function W. ( ) ( )

W W= F =ψ F (2.47)

Substituting this relation into Calusius-Duhem inequality, gives

: 0

W

∂

Σ − ≥

∂E E (2.48)

Since the Clausius-Duhem inequality holds regardless of any particular , it is necessary that 0

W

∂

Σ − =

∂E (2.49)

The second Piola-Kirchhoff stress tensor is conjugate to the Green-Lagrange strain tensor γγγγ. The constitutive equation can rewritten as

W

∂ =

∂ (2.50)

Application of the relation between the Green-Lagrange strain tensor γγγγ and the right-Green strain tensor C leads to the constitutive relation in the form

2 W W ∂ ∂ ∂ = = ∂ ∂ ∂ C C C (2.51)

Hyperelastic materials can undergo large reversible deformations. The ideal hyperelastic model of a material is based on the existence of a potential function or the strain energy function, which depends only on the terminal state of strain and in no way on the history of the straining. Thus there is a unique relation between stress and strain. The strain energy function can be expressed as a function of the principal invariants of the right Cauchy-Green strain tensor C

(

1, ,2 3

)

(23)

Chapter 2 13

The constitutive relations are obtained using the chain rule for derivation of the strain energy function. The most common expression for the constitutive equations in finite elasticity is the relation between the second Piola-Kirchhoff stress and the Cauchy-Green strain tensor, 2 i i=1, 2, 3 i I W I ∂ ∂ = ∂ ∂C (2.53) Substitution of first order derivatives of the invariants with respect to the right-Green stretch tensor C, yields 1 1 3 1 2 2 3 2 W I W W I W I I I I − ∂ ∂ ∂ ∂ = + − + ∂ ∂ I ∂ C ∂ C (2.54) We also can directly compute the Cauchy stress or true stress tensor

2 3 1 3 1 2 2 1 T 2 _I W W _I W W J J I I I I σ = = ∂ + ∂ + ∂ −∂ ∂ ∂ ∂ ∂ F F I B B (2.55) Because 2 1 1 2 3 I I I − − + − = B B I B 0 (2.56) we obtain 1 2 3 3 2 3 1 2 1 T 2 _I W _I W W _I W J J I I I I σ ₌ ₌ ∂ ₊ ∂ ₊∂ ₋ ∂ − ∂ ∂ ∂ ∂ F F I B B (2.57) Not only the mathematical invariants Ii of the tensor C can be employed as the fundamental quantity to determine the strain energy function, but also the principal stretches can be used for this purpose. The principal stretches are the eigenvalues of the deformation gradient tensor F and can be identified with the ratio of lengths of an elementary cube in the

reference configuration and its deformed representative in the actual configuration. If the strain energy function is expressed as a function of the principal stretches, the principal invariants can also be presented by means of the principal stretches λ1, λ2 and λ3

2 2 2 1 1 2 3 2 2 2 1 1 2 3 2 2 2 2 2 2 2 1 2 2 3 3 1 2 2 2 3 1 2 3 I I I I λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ − − − − = + + = + + = + + = (2.58)

The strain energy function can also be expressed as a function of the principal stretches

(

1, ,2 3

)

W W= λ λ λ (2.59)

The stress and deformation for incompressible material are related by

(24)

14 Chapter 2

where σ1, σ2, σ3 are the principal components of the Cauchy (or true) stress tensor and p is a Lagrange multiplier and is referred to as the “arbitrary” hydrostatic pressure. It arises from the constraint

1 2 3 1

λ λ λ = (2.61) This may be used to express λ3 in terms of λ1 and λ2, which are then adopted as the independent deformation measures. Accordingly, the strain energy function may be expressed as

(

1 1

)

1, ,2 1 2

W W₌ λ λ λ λ− − _(2.62)

Substitution of Eq. (2.62) into Eq. (2.60) and elimination of p yields

(

)

3 1, 2 W β β β σ σ λ β λ ∂ − = = ∂ (2.63)

For a compressible material

(25)

Chapter 3 Finite element aspects

The finite element method used for numerical calculation is described with special emphasis to the large strain analysis. If large deformations are taken into account, the equilibrium equations governing the deformation processes are nonlinear even if a linear constitutive relation is used. Analytical solutions for geometrically nonlinear problems can only be obtained for extremely simple geometries and loading conditions. Practical problems are limited to numerical solutions. A very popular numerical method for solving nonlinear mechanical problems is the finite element method.

3.1 Principle of virtual work

The basis of the displacement-based finite element method is the principle of virtual work (which we also call the principle of virtual displacement). The principle states that the equilibrium of the body requires that for any compatible small virtual displacement imposed on the body in its state of equilibrium, the total internal virtual work is equal to the total external virtual work.

int ext

W W

δ =δ (3.1)

This is the weak form of the equilibrium equation. For static mechanical problems the virtual work can be expressed as

0 0 0 0 0 0 0 T t t T T V V S dV dV dV δ +∆ ₌ ρ δ_{u g} ₊ δ_{u t} _(3.2)

where the subscript 0 refers to the undeformed configuration or the reference state, while the superscript t + ∆t indicates that the value of a quantity is considered at ‘time’ t +∆t. The superscript T denotes a transpose, the δ-symbol denotes the first variation of a quantity, ρ is the mass density and g the gravity acceleration vector, t represents the normal tractions on the

surface of the body and u is the displacement vector. The second Piola-Kirchhoff stress

tensor can be decomposed according to

t+∆t _{= + ∆}t _(3.3)

By substituting, Eq. (3.2) can be recast as

(26)

16 Chapter 3

(

)

(

)

(

) (

)

1 1 + 2 2 j j i i k k k k ij ij j i i j u u u δu δ u δu u δu γ δγ ξ ξ ξ ξ ∂ + ∂ + ∂ + ∂ + + = + ∂ ∂ ∂ ∂ (3.5)

Subtraction of γij on both sides according to the definition gives the virtual increment δγij, but now a second order term in the virtual displacement δuk emerges. Again, the incremental form of γij + ∆γij of the strain tensor can be derived. The strain increment ∆γij contains linear and quadratic contributions in the displacement increment ∆uk. The parts that are linear in the displacement increment are denoted as ∆εij and the second order contribution to the strain tensor are denoted as ∆ηij

ij ij ij γ ε η ∆ = ∆ + ∆ (3.6) with 1 ₊ 2 j i k k k k ij j i i j i j u u u u u u ε ξ ξ ξ ξ ξ ξ ∂∆ ∂∆ ∂∆ ∂ ∂ ∂∆ ∆ = + + ∂ ∂ ∂ ∂ ∂ ∂ (3.7) and 1 2 k k ij i j u u x x η ∂∆ ∂∆ ∆ = ∂ ∂ (3.8)

Introducing the virtual Lagrangian strain tensor, we have

0 0 0 0 0 0 0 0 0 0 0 0 0 0 T T t T V V V T T T t V S V dV dV dV dV dV dV δ δ δ ρ δ δ δ ∆ + + ∆ = + − u g u t (3.9)

The incremental stress-strain relation can be obtained from Eq. (2.50)

(

)

2 2 4 2 2 W W ∂ ∂ ∆ = ∆ = ∆ + ∆ = = ∂ ∂ D D D C (3.10)

The matrix D contains the instantaneous stiffness moduli of the material model. For linear

elasticity it simply reduces to Hooke’s law. Substituting this relation into Eq. (3.9) yields

(

)

(

)

0 0 0 0 0 0 0 0 0 0 0 0 0 0 T T t T V V V T T T t V S V dV dV dV dV dV dV δ δ δ ρ δ δ δ ∆ + ∆ + + ∆ + ∆ = + − D D u g u t (3.11)

This equation is linearised and consequently the second order terms in the incremental displacement are neglected. We obtain

0 0 0 0 0 0 0 0 0 0 0 0 T T t T T T t V V V S V dV dV dV dV dV δ D∆ + δ = ρ δu g + δu t − δ (3.12)

(27)

Chapter 3 17

L

∆ = ∆ =L u LH a B a ∆ = ∆ (3.13) where ∆a is the incremental nodal displacement vector, H is matrix of the interpolation polynomials for displacement or shape functions, the matrix L contains differentials which

must be determined for every strain measure and every strain condition. The BL matrix is the relation between the linear part of the stain tensor and the displacement. So, the linear contribution to the tangent stiffness in the left-hand side of the Eq. (3.12) can be expressed as

( )

0 0 0 0 T T T T L L L V V dV dV δ D∆ = δa B DB a∆ = δa K a ∆ (3.14)

Formally, the nonlinear contribution in the left-hand side of the Eq. (3.12) can be rewritten as

( )

0 0 0 0 T t T t NL NL NL V V dV dV δ = δa B T B ∆a = δa K ∆a (3.15)

In this equation the second Piola-Kirchhoff stress is represented in matrix form Tt_. Introducing those equations leads Eq. (3.12) into the form

(

KL+KNL

)

∆ =a F a (3.16) in which 0 0 0 0 T 0 T 0 0 T t 0 a L V S V dV dV dV ρ = + − F H g H t B (3.17)

3.2 Formulation for the hydrostatic stress-strain relation

For near-incompressible materials the hydrostatic pressure p is treated as a separate independent unknown to avoid overestimation of the volumetric deformation, which leads to stiffness locking since the ratio of bulk modulus over shear modulus tends to infinity. If we decompose the strain energy function into a deviatoric part Wd and hydrostatic part Wh,

2 d 2 h 2 d 2 d h W W W J p ∂ ∂ ∂ ∂ = + = + = + ∂C ∂C ∂C ∂C (3.18) in which h p W p kw J ∂ = = ∂ (3.19)

where wp is defined as a pressure function. Then the incremental stress-strain relation is obtained by differentiation of Eq. (3.18) one more time with respect to Green-Langrange strain tensor γγγγ and hydrostatic pressure p, respectively

(28)

18 Chapter 3

The right side of virtual work Eq. (3.12) should have one more term

0 0 0 0 0 0 0 0 0 0 0 0 0 0 T T t T V V V T T T t V S V dV dV pdV dV dV dV δ δ δ ρ δ δ δ ∆ + + ∆ = + − D d u g u t (3.22)

Similar to the displacement field, interpolation polynomials are introduced to describe the pressure field,

p

∆ = ∆N p (3.23)

N is the matrix of interpolation polynomials for pressure. Then, we have

0 0 0 0 T T T T L p V V pdV a dV a δ d∆ =δ B dN ∆ =p δ K p ∆ (3.24) Eq. (3.16) becomes

(

KL+KNL

)

∆ +a K p F p∆ = a (3.25)

This formula contains two sets of independent unknowns, one is the incremental nodal displacement vector ∆a and the other is the incremental pressure points vector ∆p. Therefore, an additional constraint needs to be formulated to ensure stress evolution. This constraint equation can be found from the hydrostatic stress-strain relation, which has been worked out and implemented into DIANA by Van den Bogert (1991). The hydrostatic stress-strain relation can expressed as

( )

2 ' T t t t p p J p p kw J kw C ∂ + ∆ = + ∆ ∂ (3.26) where the prime denotes differentiation with respect to J. The weak form of this relation reads

(

)

0 0 ' 0 t T t p p V p kw kw p p dV δ + d ∆ − − ∆ = (3.27)

which must hold for any admissible virtual pressure δp. Substituting the interpolation schemes for displacement and pressure fields and omitting the quadratic terms in ∆a, yields

1

T

p∆ +k− ∆ = p

K a M p F (3.28) in which the following definitions have been used:

(29)

Chapter 3 19

Bogert applies the latter method to simplify the process and save computational costs. In this study, Eq. (3.28) is solved at element level:

(

)

1 T

p p

k −

∆ =p M F K a − ∆ (3.31) Substituting this result into Eq. (3.25) yields

(

1 T

)

1

L+ NL−k p − p ∆ = a−k p − p

K K K M K a F K M F (3.32)

Once the above nonlinear system of equations is established, it is necessary to find its solution. Its character is inherently nonlinear so a nonlinear solver must be used. Basically, the nonlinear equation system can be solved by an iterative scheme. Some methods are frequently used. One of them is the Newton-Raphson method, in which two subclasses can be distinguished: the regular and the modified Newton-Raphson method. A second method is Quasi-Newton method (also called “Secant method”), where the secant stiffness matrix generally is not unique and it is illustrated in different forms, e.g., Broyden method, BFGS (Broyden-Fletcher-Goldfarb-Shanno) method and Crisfield method. Another method is called the constant stiffness method, in which the stiffness matrix will remain unchanged during this method is applied. In our coming numerical calculation, the regular Newton-Raphson method is mostly used since it has the characteristic of rapid convergence.

However, the inclusion of the pressure degrees-of –freedom on element-level may lead to spurious stress patterns in the converged displacement solutions (Sussman and Bathe (1987)). Convergence of the internal hydrostatic pressure is not necessarily reached at the same iteration step of convergence of the displacement degrees-of-freedom [∆a]. This is obvious since in the iteration process the main concern is to obtain the convergence solution of [∆a], and not the convergence of [∆p]. For more discussion about this topic see e.g. Van den Bogert (1992).

3.3 Special user subroutines for nonlinear hyperelastic materials

In order to implement new rubber models in the finite element software DIANA, a new user subroutine USRRUB has to be programmed. This routine should be capable of dealing with stress-softening analysis, so it must update

- The stress vector.

- The tangent stiffness matrix.

- The internal state variables and indicators (when we simulate material behaviour including damage phenomena).

Based on the new user subroutine USRRUB, a new interface framework for stress and stiffness are set up. This framework aims at dealing with hyper-elastic materials including or excluding stress softening phenomenon. The main steps of the framework are indicated in Fig. 3.1. This framework firstly receives general data from the main program, for example,

(30)

20 Chapter 3

coded in the Improved User-Supplied Subroutine. USRVAL are the user-supplied material parameters. If we simulate the behaviour including stress softening, we may give initial values of some internal state parameters by USRSTA. This framework finally provides the upgraded total stress vector and tangential stiffness matrix D, furthermore, the internal state variables, which may be the maximum strain energy function and the damage parameter or others if necessary, are recorded in USRSTA for later use.

A list of symbols used in Fig. 3.1.

ε0 strain vector at the start of increment

ε total strain increment USRIND user-supplied model name

USRVAL user-supplied material parameters USRSTA initial values of internal state parameters

CDM & GAO’ MODEL Continuum Damage Mechanics model combined with GAO elastic model

PEM & GAO’ MODEL Pseudo-Elastic Model combined with GAO elastic law Emax maximum strain energy

d damage variable in Continuum Damage Mechanics model dm maximum damage variable

damage variable in Pseudo-elastic model stress

0 stress for undamaged materials in Continuum Damage

Mechanics model and for loading path in Pseudo-elastic model D stiffness

D0 stiffness for undamaged materials in Continuum Damage

(31)

Chapter 3 21

(32)

(33)

Chapter 4 Constitutive models for hyper-elastic materials

4.1 Overview of constitutive models

Many attempts have been made to develop a theoretical stress-strain relation that fits experimental results for hyperelastic materials (Yeoh, 1993; Charlton and Yang, 1994; Boyce and Arruda, 2000; Miehe, 2004). There are two rather different approaches to the study of rubber elasticity. On the one hand, the phenomenological theory treats the problem from the viewpoint of continuum mechanics. This approach constructs a mathematical framework to characterize rubbery behaviour so that stress analysis and strain analysis problems may be solved without reference to microscopic structure or molecular concepts. On the other hand, the statistical or kinetic theory attempts to derive elastic properties from some idealized model of the structure of vulcanised rubber. This theory is one of the cornerstones of our understanding of the macromolecular nature of rubber.

The basic features of stress-strain behaviour have been well modelled by invariant-based and/or stretch-based continuum mechanics theories. For example, Mooney (1940) proposed a phenomenological model with two parameters based on the assumption of a linear relation between the stress and strain during simple shear deformation. Later, Treloar (1944) published a model based on the statistical theory, the so-called neo-Hookean material model with only one material parameter. However, this was proved to be merely a special case of the Mooney model. In 1950, Rivlin (1948a, b, 1949) modified the Mooney model to obtain a general expression of the strain energy function expressed in terms of strain invariants. One of the successful models in this class has recently been developed by Yeoh (1993) in the form of a third-order polynomial of the first invariant of the right Cauchy-Green tensor. An alternative high-order polynomial model of the first invariant has been proposed by Gent (1996) and takes the form of a natural logarithm.

In 1972, Ogden (1972a, b) proposed a strain energy function expressed in terms of principal stretches, which is a very general method for describing hyperelastic materials. An excellent agreement has been obtained between Ogden’s formula and Treloar’s experimental data for extensions of unfilled natural rubber up to 700% (Ogden, 1972a, b). However, the parameter identification is complicated because of the purely phenomenological character of the Ogden strain energy function.

(34)

24 Chapter 4

Gaussian network model consisting of three independent sets of chains (each set contains n/3 chains per unit volume), a so-called three-chain model. Treloar (1946) and Flory and Rehner (1943) simplified the complex polymer network by a four-chain model, in which four non-Gaussian chains are connected to the corners of a tetrahedron. Later on, Arruda and Boyce (1993) developed an eight-chain model and Wu and Van der Giessen (1993) proposed the so-called full-network model. Recently, Boyce (1996) compared the eight-chain model to the first invariant-based Gent model and demonstrated the almost equivalence of these two models in the sense of their constructions and fitting qualities of test results. A recent model, which is called the non-affine micro-sphere model has been proposed by Miehe (2004).

Theoretical analysis and engineering application require the constitutive law to be expressed as simply as possible. However, simplicity often violates rationality. When we consider a problem with a singular point (such as a crack tip or a concentrated force), the situation is different from the finite deformation case (Knowles and Sternberg, 1973; Mooney, 1940). Actually, near a singular point in rubber-like materials, the strain goes to infinity (very high value), which complicates the problem. To reflect the material behaviour near a singular point, Gao proposed a simple elastic law that separately considers the resistance of materials to tension and compression (Gao, 1997). This constitutive relation was successfully used to analyse singular problems (Gao, 2001, 1999, 1997).

In this chapter attention is restricted to the development of Gao’s constitutive model (Gao, 1997) for computational analysis of the static behaviour of hyperelastic materials. The overview of constitutive models is presented in this paragraph. In paragraph 4.2, Gao’s constitutive elastic law is presented as well as Drucker’s stability postulate and parameters are discussed; Special subroutines are programmed to implement Gao’s elastic law in the DIANA finite element package. Computational aspects are considered in the last Section.

4.1.1 Constitutive relations for incompressible materials

From the previous paragraph we have seen that, for nonlinear hyperelastic isotropic materials, the stress situation in a material point can be derived from a strain energy function. In the past decade there have been great efforts to model the mechanical behaviour of rubberlike materials and to determine the form of the elastic potential.

Numerous polymeric materials can sustain finite strains without noticeable volume changes. Since the bulk modulus of rubber is very high, several orders of magnitude larger than the shear modulus, such types of material may be regarded as incompressible so that only isochoric motions are possible. For many cases, this is a common idealization and accepted assumption often invoked in continuum and computational mechanics. For an ideally incompressible material the incompressibility constraint can be expressed in the form

det( ) 1

J = F = (4.1) The strain energy function (Eq. (2.52) can be written as

(

*

) (

*

)

1, ,2 1, 2

W W I I J= =W I I (4.2)

(35)

Chapter 4 25

When the stress-strain relation is elaborated, the assumption of an incompressible material model leads to an undefined hydrostatic pressure which can be solved by the additional constraint (4.1). The above two relations, the constitutive model and the incompressibility constraint, can be rearranged in terms of the Lagrange multiplier method yielding:

(

*

)

1, 2 ( ) W W I I= + pf J (4.3) with (1) 0 f = (4.4)

with p the Lagrange multiplier which has the same dimension as stress.

This formulation necessitates the determination of partial derivatives in the stress-strain relation following by chain differentiation. Substituting Eq. (4.3) into Eq. (2.51), the second Piola-Kirchhoff stress has the form

( )

*

( )

1, 2

=2 ∂W I I +∂pf J

∂C ∂C (4.5)

Rivlin proposed a most general form of the strain energy function in terms of the invariants for incompressible materials, namely

( )

*

(

)

(

*

)

1 2 1 2 1 , _ij 3 i 3 j i j W I I ∞ C I I + = = − − (4.6)

where Cij are material constants. Different authors consider different coefficients in this equation to derive their specific models. An overview of coefficients used in different models is given in Table 4.1.

Specifically, for incompressible materials subjected to uniaxial tension or compression, Rivlin gives the relation

(

)

(

)

(

)

* * 1 2 1 2 * 2 1 2 , ₁ , 2 t W I I W I I I I σ λ λ λ− ∂ ∂ = + ∂ ∂ − (4.7)

where σt is the engineering stress, which refers to the undeformed state. For simple shear, the relation becomes

( )

*

( )

* 1 2 1 2 * 1 2 , , 2 s s W I I W I I I I τ _µ γ ∂ ∂ = + = ∂ ∂ (4.8)

where τs is the shear stress and γs is the shear strain, which is related to the invariant I1 by

γs2=( I1-3), µ is the shear modulus.

In the above-mentioned models, two very important constitutive models should be emphasized. One is the neo-Hookean model and the other is the Mooney-Rivlin model. Taking only the first term of Eq. (4.5) it yields the neo-Hookean model

(

)

10 1 3

(36)

26 Chapter 4

Table 4.1 The coefficients used in different material models based on the strain energy function expressed as an infinite series in terms of I1 and I2*.

Model reference coefficients Neo-Hookean (Treloar, 1944) C10

Mooney-Rivlin (Mooney, 1940) C10 C01

Tschoegl (Tschoegl, 1971) C10 C01 C11

Isihara et al. (Isihara et al., 1951) C10 C01 C20

James et al. (James et al., 1975) C10 C01 C11 C20 C02

Tschoegl (Tschoegl, 1971) C10 C01 C22

Yeoh (Yeoh, 1993) C10 C20 C30

Biderman (Traloar, 1975) C10 C01 C20 C30

James et al. (James et al., 1975) C10 C01 C11 C20 C30

Haupt and Sedlan

(Haupt and Sedlan, 2001) C10 C01 C11 C02 C30

James et al. (James et al., 1975) C10 C01 C11 C20 C02 C21 C12 C30 C03

James et al. (James et al., 1975) C10 C01 C11 C20 C02 C21 C30 C03 C40

Lion (Lion, 1997) C10 C01 C50

Treloar (1944) constructed the strain energy function on the bases of Gaussian statistics and a molecular network theory

(

1

)

1

3 2

W = NkT I − (4.10)

where N is the number of network chains per unit volume, k is Boltzmann’s constant and T is the absolute temperature.

Note that Eq. (4.10) is the Gaussian statistical model, which is equivalent to the continuum mechanics model Eq. (4.9), even though the phenomenological and statistical theories start from completely different premises. A shear modulus µ =C10 is the only

parameter of the neo-Hookean model for the incompressible material. This simple form provides an adequate first approximation to the behaviour of rubber-like solids and is generally regarded as a valid prototype for this class of materials. It has the advantage of being easy to treat mathematically so that solutions to many problems have been obtained by its use.

(37)

Chapter 4 27

(

2 2 2

)

(

)

(

*

)

10 1 2 3 01 2 2 2 10 1 01 2 1 2 3 1 1 1 3 3 3 3 W C λ λ λ C C I C I λ λ λ = + + − + + + − = − + − (4.11)

When simple shear case is considered, we have

(

10 01

)

* 1 2 2 2 s s W W C C I I τ _µ γ ∂ ∂ = + = + = ∂ ∂ (4.12)

The Mooney-Rivlin and neo-Hookean strain energy functions have played an important role in the development of the non-linear hyperelastic theory and its applications (Ogden, 2001, 1984). It has been proved by Bogert (1991) that the Mooney model performs well for moderately large deformation of uniaxial elongation and shear deformation. But, it cannot describe the S-curvature of the force-stretch relation in the uniaxial elongation experiment and the force-shear displacement relation in a shear experiment. Tschooegl (1971) suggests that failure of the Mooney-Rivlin model to provide adequate multiaxial data predictions arises from not including enough terms of the possible expansions of Eq. (4.6). Rivlin and Saunders (1951) suggested that C01 should not be constant, but should be dependent on I2,

decreasing with increasing I2.

Working within the continuum mechanics framework for the strain energy function as proposed by Rivlin, Eq. (4.6), several investigators have used higher order terms in I1 and, in

some case I2, to account for the departure from neo-Hookean/Gaussian behaviour at large

stretches. One of this type is the Yeoh’s model (1993) according to

(

)

(

)

2

(

)

3

10 1 3 20 1 3 30 1 3

W C I= − +C I − +C I − (4.13)

Comparison of Yeoh’s model with the published biaxial data of James and Green (1975) shows a good ability to predict multiaxial data. Using the higher order terms I1 in the strain

energy function has been shown to work well in capturing different deformation states at moderate to large deformations. An alternative high order I1 model has recently been

proposed by Gent (1996) and takes the form:

1 3 ln 1 6 _m 3 I E W I − = − − − (4.14)

where E is the small strain tensile modulus and the Im is the maximum value for I1 where as I1

approaches Im the material approaches limiting extensibility. As discussed in Boyce (2000), the natural logarithm term in the Gent model can be expanded to yield the following expression for the strain energy:

(

_{) ( )( ) ( ) ( )}

(

)(

)

(

)

2 3 1 1 2 1 1 1 1 1 3 3 3 6 2 3 ₃ ₃ 1 3 1 3 m m n n m E W I I I I _I I n I + = − − + − + − + ⋅⋅⋅ − − + − + − (4.15)

(38)

28 Chapter 4

Apart from stretch invariants also the principal stretches can be used as the set of basic parameters that describe the material behaviour for incompressible material. Varga (1966) proposed the first strain energy function as the sum of the principal stretches.

(

1 2 3 3

)

V

W =µ λ λ λ+ + − (4.16)

in which µv is a material parameter.

He has demonstrated that the range of validity of this strain energy function is comparable with the neo-Hookean strain energy function. It should be reminded that this model could not be expressed as a simple function of the invariants I1 and I2*. More general, Valanis and

Landel (1967) proposed that the strain energy function for an incompressible material could be expressed as a symmetrical and separable function of principal stretches

( ) ( ) ( )

1 2 3

W w= λ +w λ +w λ (4.17)

To be able to predict an accurate stress state at really large strains Ogden (1972a) proposed a strain energy function in terms of the principal stretches for incompressible materials

(

1 2 3

)

1 3 r r r r r r W µ λ λ λα α α α ∞ = = + + − (4.18)

in which µr and αr are material parameters. The indices αr need not to be integers.

Applying the general relation for principal Cauchy stresses to the Ogden material model we obtain

(

)

1 1, 2,3; no sum r i r i r p i α σ ∞ µ λ = = − = (4.19)

where p is an arbitrary hydrostatic pressure introduced because of the incompressibility constraint.

We can make some more general statements about the nature of the coefficients. Consideration of stability and a physically realistic response lead to the constraint (see paper (Ogden 1972a, Yeoh, 1997)) to every pair of coefficients

0

r r

µ α > (4.20)

Ogden also mentioned that two terms are sufficient to describe the load-stretch characteristics in tension and shear, but the third term is necessary to include the behaviour in biaxial deformation. The value of α1 is usually less than 2 (Yeoh 1996). This reflects the

initial decrease in shear modulus with increasing strain. The value of α2 is greater than 2 as it

is necessary to account for the subsequent increase in shear modulus at large strains.

4.1.2 Constitutive relations for compressible materials

(39)

Chapter 4 29

hand various kinds of biological tissues and synthetic elastomers are hyperelastic and significantly compressible. Besides the mentioned incompressible form of W, which, in fact, only describes the distortional behaviour of rubber-like materials, additional terms in the strain energy function can be introduced to deal with pure volumetric deformations. For such a model a modified set of matrix invariants i, which are a function of the regular I-invariants,

are applied to achieve a strict separation of the strain energy function in a deviatoric part and a volumetric part:

(

1 2

)

( )

3

d v

W W I= +I +W I (4.21)

where Wd is the distortional part of the strain energy function, described in the preceding paragraph and Wv is the volumetric part of the strain energy function. The distortional (volume-preserving) deformation is described by a modified deformation gradient

1/ 3

J−

=

F F (4.22)

and we define the modified right Cauchy-Green strain tensor as

2/3

J−

=

C C (4.23)

The regular invariants of right Cauchy-Green strain tensor C, denoted as Ii, should be identical to the invariants of J2/3C , in which the modified invariants of C will be denoted as

i. We have

(

)

2 2/ 3 4/ 3

1 1 2 2 and 3 det 1

I ₌J− I I ₌J− I I ₌ _F ₌ _(4.24)

This approach does not predict the volume change with sufficient accuracy in uniaxial elongation experiments. Van den Bogert (1992) introduced a coupling term Wc to improve the results. Wc depends on distortional as well as on volumetric deformations

d v c

W W= +W W+ (4.25)

Computational modelling of rubber-like materials under monotonic and cyclic loading